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1. INTRODUCTION 


1.1 Background of the Present Effort 

Design optimization is an issue of continuing importance in a variety of 
aeronautical applications, including rotorcraft performance analysis. Developing an 
effective compromise design that can meet the requirements of rotor performance in 
hover and in climb or high speed axial flight is an exceptionally difficult problem. Hover 
itself remains an important design point because of the critical role it plays in sizing the 
engine and establishing rotor disk loading and the blade twist distribution. Moreover, the 
axial flight condition has become even more critical, since the success of present and 
proposed tiltrotor designs depends on efficient performance of lifting rotors in propeller 
mode. The complexity of the problem increases still further when structural and dynamic 
constraints are considered. The primary objective of this effort was to develop a 
comprehensive design optimization capability suitable for analyzing both hover and axial 
flight conditions in modem helicopter and tiltrotor configurations. 

An essential prerequisite for a practical and effective optimization analysis is for a 
refined, flexible, and well-validated hover performance code to be used as the basis of the 
performance calculation. The task of developing computational tools to reliably analyze 
rotor performance has been under study for many years. Much of the work in this area 
has centered around efforts to approximate crucial rotor wake effects through 
generalizations of vortex wake trajectories from empirical data (Refs. 1 and 2) or through 
interpolations between 'representative' free vortex computations (Ref. 3). It has long been 
clear that the preferred approach is the explicit computation of full free wake vortex 
flows; this approach avoids reliance on particular data sets and yields force-free, 
physically valid wakes which enhance confidence in performance predictions. The next 
section will summarize previous work in the development of free wake analyses as well 
as describe the formulation implemented here. 

A second requirement for a successful optimization analysis is that it must 
efficiently evaluate candidate designs with minimal resort to repeated calls of the 
computationally expensive open-loop performance analysis. It is in this regard that the 
influence coefficient approach proves highly advantageous. An analysis incorporating 
this approach, designated EHPIC (Evaluation of Hover Performance using Influence 
Coefficients), has proved to be the ideal foundation for a free wake optimization routine, 
since its fomiulation permits straightforward evaluation and exploitation of information 
on the gradients in performance due to design changes. The results presented below will 
show how an optimization analysis based on influence coefficients can perform efficient 
exploration of design space without sacrificing the refined physical model associated with 
free wake treatments. 

Previous reports and papers (Refs. 4 and 5) have described exploratory work on 
this topic, involving computations of improved performance on representative rotor 
configurations using a loose coupling of the EHPIC code to a linear optimization 
analysis. While these demonstration calculations were encouraging evidence of the 
feasibility of developing an optimization code based on EHPIC, they were restricted to 
using twist as a design variable and lacked the flexibility to realistically constrain the 
evolving designs. It was thus recognized that considerable additional development and 
validation work was required to implement a useful, practical optimization code. 
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References 6-9 summarize extensive performance validation efforts undertaken 
with the EHPIC code. These references describe correlation studies involving over thirty 
separate rotor configurations, including tiltrotor designs, conventional main rotors, and 
tail rotors. In general, these studies achieved quite accurate correlation with minimal 
resort to "dialing" of model parameters. This work in OGE (Out of Ground Effect) hover 
calculations has been supplemented by enhancement of the original EHPIC (Mod 0.0) 
code (Ref. 9) to include a ground plane to allow IGE (In Ground Effect) calculations 
(Ref. 10). Validation of IGE performance and wake geometry predictions discussed in 
References 10 and 11 has also been encouraging. 

Thus, the underlying EHPIC hover model has undergone substantial testing and 
refinement since the original work described in Reference 9, but still other improvements 
were judged desirable to improve the flexibility and physical fidelity of the code. Many 
of these new features are described in Sections 2, with additional details on 
implementation and use provided in Reference 12. The new features of the model 
include: a provision for lift limitation due to airfoil stall; direct calculation of wake- 
induced velocity at specified points; refinements in vortex core modeling; research on 
vortex/blade encounter loading; and significant CPU reduction through the streamlining 
of influence coefficient evaluations. Sections 3 and 4 describe still more substantial 
improvements to the baseline EHPIC analysis, including a provision for high-resolution 
roll-up of tip vortices and computation of static structural deflection. Later sections 
direcdy take up the topic of the optimization algorithms implemented as well as the 
calculations undertaken to validate the functionality and capabilities of these algorithms. 

1.2 Review of Previous Work in Rotor Design Optimization 

There exists a large body of literature on numerical optimization procedures. 
Recent reviews of aerospace- and rotorcraft-oriented applications of numerical 
optimization give evidence of the maturity of these techniques (Refs. 13-14). As for 
particular applications, Moffitt and Bissell (Ref. 15) undertook a general examination of 
rotor airloads in both hover and forward flight using a prescribed wake aerodynamic 
model with circulation coupling. Their aim was to find airload distributions that led to a 
minimum power required for a specified thrust Nagashima and Nakanishi (Ref. 16) 
studied hover performance optimization for coaxial rotors, employing both a closed-form 
"generalized momentum" model of the wake as well as a simplified free wake model 
based on vortex rings. Walsh, Bingham, and Riley (Ref. 17) and Chattapodhyay, Walsh 
and Riley (Ref. 18) describe the assembly of several existing aerodynamic and vehicle 
trim models and a commonly available optimization routine (CONMIN) into a broadly 
applicable blade design optimizer for hover and forward flight This work is part of a 
larger effort to develop a multidisciplinary design optimization for rotorcraft (Adelman 
and Mantay, Ref. 19) which has been directed at achieving appropriate compromise 
designs across a wide range of performance measures. This in turn has been 
supplemented by more specialized studies performed by Chattapodhyay and McCarthy 
(Ref. 20). 

Nearly all of the studies mentioned above rely on qualitatively similar approaches; 
an objective function (e.g., vibratory load levels or hover figure of merit) is defined along 
with a set of design parameters to be varied. Limitations are placed on the variations 
allowed in these quantities in the form of inequality constraints and initial estimates of 
parameter values to meet the specified performance levels are made. Using 
approximations to the actual nonlinear relations governing the objective function, the 
design parameters are varied in the vicinity of the initial estimate to reach the desired 
level of this function. The more recent work, (Refs. 19 and 20), has expanded the scope 
of such methods to address multidisciplinary optimization efforts involving the 
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simultaneous maximization or minimization of complex objective functions composed of 
weighted measures of performance (e.g. hover figure of merit and blade dynamic 
frequency placement, or climb power and noise signature). 

All this work represents a broad range of accomplishment in design optimization, 
although as a general rule these efforts have featured simplified aerodynamic models in 
their performance calculations (e.g., strip theory and prescribed wake models). This 
doubtless has been a response to the need to provide designers with computationally 
efficient tools; however, now that more and more computational power is becoming 
routinely available, it is appropriate to employ more advanced methods. Simplified 
models lack the generality, accuracy and adaptability of the more advanced treatments 
now available and so their use in an optimization analysis undermines confidence in the 
true optimality of the designs computed. The advantage of the cuirent effort is that the 
advanced free wake analysis embodied in the EHPIC code can be utilized to produce a 
more comprehensive, reliable and flexible optimization treatment Furthermore, EHPIC's 
influence coefficient approach helps ameliorate potential increases in CPU tim e because 
in the course of its own solution it calculates many of the derivatives required by the 
optimization algorithm. 
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2. HOVER PERFORMANCE PREDICTIONS USING A FREE WAKE ANALYSIS 


This section recapitulates the development of the EHPIC code and introduces the 
important concepts that were carried over into the cunent effort. The evolution of the 
EHPIC code from the Mod 0.0 version originally delivered to NASA in 1987 to the later 
Mod 1.0 and Mod 2.0 versions is summarized. Also described in this section are the 
major modifications that were incorporated in EHPIC/HERO to enhance performance and 
to adapt the EHPIC analysis to the tasks associated with design optimization. 

2. 1 The Influence Coefficient Approach to Free Wake Analysis 

Early efforts to develop free wake hover models using time domain calculations 
(Refs. 21 and 22) were hampered by long computation time and poor (or nonexistent) 
convergence due to the inherent instability of the hovering rotor wake (Fig. 2-1). 
Reference 9 describes the development of an influence coefficient relaxation approach to 
the free wake problem that determines the wake geometry while circumventing the 
convergence problems associated with the time-marching simulations. As mentioned 
previously, the EHPIC code has produced accurate performance predictions for a wide 
variety of rotor systems in hover, and has proved to be flexible and robust. EHPIC also 
includes curved vortex elements (Reference 23) in the model of the rotor wake to enhance 
both the efficiency and the accuracy of the computation. 

The general objective of a free wake hover analysis is to find the wake geometry 
that satisfies two conditions: first, that the wake filaments are in free motion; and second, 
that the flow tangency condition is satisfied on the blade. To achieve the free motion 
condition, the wake filament trajectories must be tangent to the local velocity vector 
evaluated on the filament when viewed in a rotating reference frame, i.e., there will be no 
crossflow velocity components at any point on the filaments under force free conditions. . 
The coupled free wake/lifting surface hover analysis in EHPIC proceeds by first making 
an initial guess for the blade loads and the wake geometry. This initial guess will not, in 
general, satisfy the required conditions, and so must be adjusted in a succession of 
solution steps. To accomplish this, the independent variables in the problem (the bound 
circulation at stations along the blade and the vortex wake position coordinates) are 
systematically perturbed, and the effect of these perturbations on the dependent variables 
(the downwash on the blade and crossflow velocities in the wake) are summed and 
formed into influence coefficients (Fig. 2-2). These coefficients allow the construction of 
a set of simultaneous linear equations in matrix form which predict the change in 
dependent variables due to the changes in independent variables. The coefficient array so 
formed can be used to null the crossflow and downwash velocities by inverting it and 
multiplying it by the vector of residual velocities. The influence coefficient array appears 
in a linear system of equations in the following form: 


Aq 


QqxQqy 


Ax 

Aw 


QwxQwy 


Ay 


( 2 - 1 ) 


The independent variables on the right hand side are, respectively, the position 
perturbations and bound circulation perturbations from the initial state, while the 
dependent variables on the left hand side represent the crossflow velocities in the wake 
and the downwash at the blade. 
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Figure 2-1. Typical form of time domain instabilities observed in free vortex 
calculations of hovering rotor wakes (dashed line); solid line 
represents the idealized contracting wake solution for this one-bladed 
rotor. 





Figure 2-2, Method of displacing a vortex to eliminate velocity components in 
planes normal to the vortex curve. Equilibrium state is zero 
crossflow velocity in every plane. 


5 



These variables, of course, represent perturbations about an initial guess, which in 
general features some nonzero residual velocities on the left hand side. Were the problem 
purely linear, inverting this matrix and multiplying it by the residual velocity vector 
would yield a vector of wake displacements and circulation perturbations that would 
exactly null the velocities in question, i.e.. 


r Ax 


Qqx 

Qqy 

-1 

flics 

i — 

3 4 
1 


Qwx 

Qwy _ 


. w rcs „ 


In general, the process must be repeated due to the inherent nonlinearity of the problem, 
and only small fractions of the residual velocities are nulled in each iteration. In practice, 
this approach has been found to have robust convergence properties, despite the 
complexity of the relaxation procedure. For most rotor configurations with up to ten 
wake filaments per blade, convergence is achieved after ten to twenty relaxation steps 
(Fig. 2-3). These solutions, while they are self-preserving vortex flows, can be shown to 
possess time-domain instabilities that preclude the success of time-marching calculations 
(Ref. 9). 

In the EHPIC code, the wake vortex filaments are represented by curved vortex 
elements which provide an efficient means to perform the Biot-Savart integrations 
necessary for the evaluation of wake-induced velocities (Fig. 2-4). A vortex lattice/lifting 
surface analysis is used to evaluate the thrust and induced drag on the rotor. An array of 
thin lifting vortex quadrilaterals are used to model each blade, and the relaxation solution 
produces a vector of bound circulation values that null the downwash at the control point 
at the center of each quad. Using the local values of free stream and induced velocity, the 
Joukowski law is applied to find the force and moment on each lattice element. This 
procedure also produces a spanwise lift coefficient distribution that is used, along with 
the local Mach number, in a look-up scheme that provides the profile drag coefficient of 
that section. 

One advantage of the influence coefficient approach is that it finds the physically 
correct, self-preserving wake geometry without the instabilities and consequent lack of 
convergence of earlier methods that convected the wake in a time-marching manner. 
Furthermore, once a converged solution has been obtained, adjacent solutions that are 
almost linearly close along a performance curve are readily obtained with only a few 
relaxation steps, eliminating the need to wash out the transients occurring in time- 
marching schemes. 

2.2 Evolution of the EHPIC Code 

The EHPIC code has evolved considerably since the original work summarized in 
Reference 9. This initial effort was directly purely at obtaining performance solutions for 
isolated rotors operating out of ground effect, and the resulting code was designated the 
Mod 0.0 version. Subsequent development work sponsored by NASA to upgrade the 
original code led to three major modifications: implementation of a ground plane model 
to allow for computations in ground effect; incorporation of an eigenanalysis package to 
permit evaluation of the linearized time domain stability of converged configurations; and 
substantial revision of the original coding to reduce the code's CPU requirements on 
vector processing computers. This version was designated Mod 1.0, and its development 
is described in detail in Reference 10. 
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Figure 2-3. Typical converged rotor wake configuration for the EHPIC code (six 
free wake filaments trailing from each blade). (Wake from only one 
blade shown.) 





POINT OF / 

EVALUATION 

Figure 2-4. Typical arrangement of elements used to discretize trailing vortex 
filaments for wake-on-wake velocity calculations. The formulation 
of parabolic Basic Curved Vortex Elements (BCVE's) and Self- 
Induction Vortex Elements (SIVE's) given in Reference 23. 
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The ground plane model implemented in Mod 1.0 consists of an image system of 
vortices located symmetrically opposite the free vortex elements in the physical wake. A 
prescribed wake model completes the transition from the freely distorting filaments to an 
efficient far wake representation. References 10 and 11 both contain discussions of wake 
geometry and performance results obtained with the ground effect model. 

This vectorized EHPIC Mod 1.0 variant achieved roughly a factor of four in 
computation time reduction over Mod 0.0 when implemented on a CRAY Y-MP. 
Additional computation time reductions - independent of computational platform - have 
been realized in a Mod 2.0 version developed under an internal effort at Continuum 
Dynamics. This version incorporates a broad range of efficiency enhancements, 
primarily limiting the frequency of updates of the influence coefficient array. Mod 2.0 
contains the same physical model as Mod 1.0, but can run three to four times faster on the 
same machine. The speed-up in Mod 2.0 is independent of (serial or vector) platform, so 
a Mod 2.0 calculation operating on a CRAY Y-MP or similar vector processor can run 
over an order of magnitude faster than the original EHPIC Mod 0.0. 

The Mod 2.0 variant was the starting point for the development of EHPIC/HERO 
described here. (Unless otherwise specified, any subsequent references to 'the EHPIC 
code’ should be understood to refer to Mod 2.0). During the evolution of the 
EHPIC/HERO code from the baseline EHPIC analysis, a variety of extensions and 
revisions were incorporated independent of any optimization functions to relieve earlier 
limitations on the analysis and effect improvements in the consistency of the model, or 
simply to make the analysis more convenient to use. The remainder of this section 
describes the implementation of the most significant of these modifications, while 
discussion of two major extensions - the calculation of tip vortex roll-up and the 
implementation of a structural deformation model - is reserved for subsequent sections. 


2.3 Calculation of Aerodynamic Loading 
2.3.1 Vortex lattice model 

In EHPIC and EHPIC/HERO, thrust and induced torque are computed using a 
vortex lattice formulation similar to that described in Reference 24. The present model 
allows substantial flexibility in the specification of the blade's planform so that complex 
designs may be accommodated. Currently, the lattice can be divided into as many as 
fifteen different regions, with separate linear distributions of twist, taper, and sweep 
within each. The spacing of the quadrilaterals in a vortex lattice analysis is an important 
consideration, as is discussed in Reference 25. The judicious selection of the density, 
spacing, and orientation of the quadrilaterals can considerably enhance the efficiency and 
rate of convergence of the blade loading. The current analysis has been provided with 
sufficient flexibility to arrange essentially arbitrary chordwise and spanwise distributions 
of lattice elements though the control points are always assumed to lie at the geometric 
center of the quadrilateral. All of the calculations discussed in this report feature uniform 
spacing of vortex quadrilaterals both in the chordwise and spanwise directions, unless 
otherwise noted. (Note: A special semicircle cosine spacing option similar to the Lan 
Type-B model discussed in Ref. 25 is currently available within the EHPIC/HERO code. 
The description of the vortex lattice structure in this section does not apply to this option. 
The semicircle spacing option was not used for any of the results presented in later 
sections.) 

The vortex quadrilateral lattice is drawn in blade coordinates, which have their 
origin at the rotor hub with Z down along the shaft, X radially outward along the 
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planform, and Y normal to the XZ plane (Figure 2-5). First the blade segments are laid 
out separately in the XY-plane applying taper and sweep. The lattice is displaced toward 
the trailing edge by a distance of one quarter of the chordwise length of the leading edge 
quadrilaterals. For one row of quads and an unswept rectangular planform, this puts the 
leading edge quadrilateral along the quarter chord line of the blade and the vortex lattice 
control points (center points of each quad) along the 3/4 chord line of the blade. The 
lattice is inset from the blade root and tip by a distance equal to a quarter of the width of 
the last quad at either edge (Ref. 26). For reference purposes, the quarter-chord line of 
the blade is taken as the tine that connects the quarter chord points of each blade section, 
while the X axis of the blade coordinate frame is the line connecting the hub with the 
quarter-chord of the root section. The sweep angle for any segment is defined as the 
angle the local quarter-chord line makes with the X axis. Pitching moment calculations 
use the local quarter-chord line as a reference axis, however collective pitch is applied 
about the X axis. 

In actual calculations, the order of operations is different than is shown in Figure 
2-5; first, taper is applied linearly from root to tip along each segment Then sweep is 
applied by displacing each segment toward its trailing edge (+Y-direction); the sweep 
angle is the angle between the X-axis and the quarter chord line. The twist gradient is 
applied by rotating each chord of the lattice about its quarter chord point. Finally, 
anhedral is applied by displacing each segment downward in the +Z-direction. The 
resulting vortex lattice structure is stored and written to a file that can be used in a 
graphical verification of the planform. Also, once assembly of the lattice is complete, it 
is subjected to a local stretching to account for the effects of compressibility on sectional 
lift and moment properties, as described in Reference 9; compressibility effects on 
sectional drag are captured through the look-up procedure described in Section 2.3.2. 

In addition to the geometric inputs just described, a camber distribution for the 
blade may also be specified. Numerical input is used to describe the geometry of the 
camber line. When camber is present, the lattice itself is not deformed to fit the specified 
distribution, rather the boundary conditions at the vortex quadrilateral control points are 
altered to introduce the surface slope into the calculation. However, a certain minimum 
chordwise density of quadrilaterals is required to resolve the camber distribution; an 
absolute minimum of three quadrilaterals chordwise should be used, with five or more 
being desirable. This level of quadrilateral density can create a substantial computational 
burden. An alternative approach for including empirical zero lift angles is also available 
for cases where improved computational efficiency is desired. An appropriate rotation of 
the vector normal to the blade at the control points of each quadrilateral will alter the 
effective angle of attack of the section and can be used to introduce the shift of the zero 
lift angle of attack. In this manner, realistic zero lift angles of attack can be introduced 
into the calculations with only one quadrilateral chordwise. The introduction of the zero 
lift angle and the option for invoking the pitching moment coefficient from two- 
dimensional tables are described later. 

The solution method used to find the bound circulation given the vortex lattice is 
essentially a straightforward implementation of the classical approach described in the 
literature on lattice methods for fixed wing applications (e.g.. Ref. 24). Each of the 
quadrilaterals is examined individually and a mean vector normal to the quadrilateral 
surface is established as shown in Figure 2-6, which also shows the location of the 
'control point' associated with the quadrilateral. Given this and the location and 
orientation of each of the quadrilaterals on the blade, the velocity induced by the blade 
lattice on each of the control points is determined, assuming unit strength for each 
quadrilateral. Then the resulting velocity is resolved in the normal direction at each 


9 


Top 




Side 


Front 


Side 


AY 


X 


a) Baseline 



b) Baseline + twist 


c) Baseline + taper 



d) Baseline + sweep 


Front 


Figure 2-5. Sample blade layout: illustration of definition of taper, twist and 
sweep on a 5 x 20 lattice. 
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Normal vector at X c is n = Rj x R 2 
X c = 0.25 (x 1 + X 2 + X 3 + X 4 ) 


Figure 2-6. Typical vortex quadrilateral, showing the comer indices and 
diagonal vectors. Control point found using the mean of the 
comers. Normal vector defined as indicated. 



a) implementation in EHPIC 



b) new implementation in EHPIC/HERO 


Figure 2-7. Schematic of the blade/wake junction with a single trailing filament. 
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control point, yielding an array of influence coefficients relating the vector of bound 
circulations y to the vector of downwashes jy at each control point: 


w = Ay (2-3) 

The vector y can then be used to solve for the forces and moments on each vortex 
quadrilateral edge by applying the Joukowski Law , i.e.. 


Ejk = P Yj Gfck x Sjk ) ljk j = 1. -n k=l,2,3,4 (2-4) 

Here, Sjk is the unit vector directed along edge k of quadrilateral j ; lj^ is the 

length of this side, and Yj is the strength of the quad. The reference velocity q for the 

evaluation of forces is computed at the midpoint of the edge k. The forces on each quad 
edge are then summed to yield the integrated forces on each blade: 

E = (HI+YI+TK) = X X Ejk ( 2 ' 5 ) 

j * 

Moments exerted by the blade about the sectional quarter-chord reference axis can also 
be computed: 

M = (LI + MI+ NS) = X X Ejk x Ijk (2-6) 

j k 


Here, i is the vector from the reference axis to the point of action of E- Moments about 
the blade root are taken to compute the coning and in-plane deflection, while pitching 
moments can be computed for each section to use for the calculation of torsional 
deformation. 

It is appropriate here to note a change implemented in EHPICVHERO involving 
the blade/wake junction. In EHPIC, the trailing vortex filaments were assumed to be 
attached to the control points at the trailing edge of the lattice. As described in the 
outline of the near wake model in Reference 9, the influence of the wake elements 
attached to the blade are deleted and replaced with an 'overlap' region consisting of an 
extension of the bound vortex lattice into the near wake. One new aspect of the 
aerodynamic model in EHPIC/HERO was the implementation of a refined model of the 
near wake that could be used to replace the overlap model. One step required for the 
implementation of this high-resolution near wake model (described in the next section) 
was to arrange a new blade/wake junction, illustrated in Figure 2-7. The filaments now 
attach to the trailing edge of the vortex lattice; this shift in the release points of the wake 
Can have a noticeable effect on rotor performance, as will be shown in Section 6. The 
overlap model used in EHPIC may still be invoked, if desired. 

2.3.2 Sectional drag and moment characteristics 

To introduce forces generated by profile (viscous and pressure) drag into the 
calculation, the only practical approach at present is the use of two-dimensional airfoil 
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data. In the EHPIC code, the form of the data input required involved tabulated values 
for sectional drag coefficients as a function of Mach number and lift coefficient. 
EHPIC/HERO now requires C-81 format tables of lift, drag, and moment coefficients as a 
function of angle of attack. However, section lift coefficient and Mach number are still 
the variables used to enter the tables, while angle of attack is simply an intermediate 
parameter. In the present calculation, the section lift coefficient is computed using the 
vortex lattice analysis described above, and the C-81 lift coefficient table is entered with 
this value and the Mach number to find an effective angle of attack for the section. Using 
this angle as a parameter, the drag coefficient (and, if desired, the moment coefficient) are 
computed for the local Mach number. The lift coefficient table is also used to find the 
zero lift angle as well as the maximum section c\ for use in the lift stall calculation 
described in Section 2.4. 

The two-dimensional coefficients are, of course, defined only for a specific airfoil 
section. Many rotor blades feature more than one section along the span. In the current 
analysis, as many as ten different sections along the span may be specified. For any 
given spanwise station, the section coefficients are computed for each of the two airfoils 
that bound the segment containing the station of interest and then interpolated linearly to 
the desired point. Similarly, for each Mach number/lift coefficient pair, bilinear 
interpolation is used to find the appropriate coefficients within each look-up table of drag 
and moment coefficients. The user has the option of applying pitching moment 
computations from either the look-up process or through the vortex lattice calculation. 

2.4 Modeling of Lift Stall 

The original EHPIC model included no provision for limiting the lift generated by 
the vortex lattice aerodynamic model. EHPIC/HERO incorporates an option to use the 
empirical information on ci max to preclude unrealistically high lift values. The present 
stall model works by tilting the vectors normal to the blade surface of quadrilaterals in 
stalled sections by the amount required to reduce ci computed without stall to the ci m „ for 

the section. This "stall adjustment angle", Aoc,, is assigned to each blade quadrilateral of 
the stalled section. Once a section is stalled, the correction angle adjusts to ensure that c\ 
never exceeds ci m „. As the calculation proceeds, it is possible for the section to drop out 
of stall, especially during the design optimization process. In this case, the adjustment 
angle resets to zero and the section is no longer considered to be in stall. 

Even though ci never increases above cj max , it is appropriate for the profile drag to 

continue to increase as the angle of attack of a stalled section increases. The value of Acts 
is exactly the increment of angle of attack of the lattice section over the stall angle. 

Therefore, when invoking the look-up tables, an angle of attack of a + A a* is used to 
compute the drag coefficient One simplification in the present model is that the section 
lift coefficient is assumed to stay constant at C| mmx as the angle of attack is increased 
beyond stall, rather than decreasing as is typically the case for 2D airfoils. The existence 
of a multi-valued lift function was found to cause dithering in the optimization 
calculations for very heavily loaded rotors. Results of demonstration calculations of the 
stall feature on a realistic rotor configuration is shown in Section 6. 

2.5 Scan Plane Calculations 

EHPIC/HERO includes a provision for computing the wake-induced velocity at 
specified sets of scan points in the vicinity of the rotor. The present implementation sets 
up points lying in "scan planes" at particular azimuth angles relative to the blade. Radial 
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and vertical spacing of points are defined by the user. Thus, planes of "crossflow" 
velocity data at a given azimuthal location relative to the blade are easily captured, while 
planes of data normal to the rotor disk can be captured simply by using results from 
within a large number of such crossflow planes at a particular vertical location. It should 
be noted that the wake geometry and thus the resulting flow field are at all times assumed 
to be axisymmetric. 


Sample calculations were undertaken to demonstrate and validate the scan plane 
velocity computations. The data set selected involved the downwash velocity 
measurements on the ,658-scalc V-22 rotor described in Ref. 27. For this case, 
EHPIC/HERO was run with nine filaments trailing from the span. Downwash 
computations were made in four crossflow planes downstream of the reference blade and 
then suitably averaged in an attempt to duplicate the time-averaged measurements given 
in Ref. 27. Two operating conditions were considered, both with a tip speed of 232 m/s 
(760 fps) and thrust coefficients 0.0117 and .0059. The comparisons of the predictions 
with the time-averaged data are shown in Figures 2-8 and 2-9. As is evident, the 
predictions capture the behavior of the measured wake reasonably well, with the 
exception of a single anomalous point in the center of the measured distribution. 

2.6 Determination of Trailing Vortex Strength 

Figures 2-10 and 2-11 illustrate the procedure for determining vortex filament 
strengths in the EHPIC/HERO code. The circulation distribution along the blade is 
divided up into zones which each correspond to a particular vortex filament The change 
in circulation across each zone is assigned to the vortex filament associated with that 
zone, as shown. In order to have a physically consistent wake model, it is necessary that 
each filament trail from the centroid of the circulation distribution it represents. 

In the original EHPIC code, the zone boundaries were input by the user and 
remained fixed throughout the calculation. After a calculation was complete, the user 
typically had to adjust the zone boundaries and vortex release points to ensure that all the 
circulation generated on the blade was trailed into the wake and that vortex filaments 
were released from the appropriate centroid positions. Though this scheme works fine for 
performance calculations at particular operating conditions, it was found to be inadequate 
for the optimization calculations to be undertaken in this effort During the optimization 
process, substantial design changes inevitably lead to changes in the circulation 
distribution that require the position and strength of the trailing filaments to adjust. 
Without some form of internal adjustment, the optimization process will almost 
invariably lead to a non-physical solution. 

The EHPIC/HERO code offers two new approaches that allow the wake model to 
adapt to a varying circulation distribution. These methods are illustrated in Figures 2-10 
and 2-11. In both methods, the code determines the circulation zones internally, updating 
zone boundaries when necessary as the circulation distribution evolves. In the first 
method, the zone boundaries are placed equidistant between vortex filament release 
points. Zones on either side of the peak circulation location are treated separately. 
Hence, if the peak circulation location moves, the zone boundaries will adjust even if the 
vortex filament release points are fixed, (as shown in Figure 2- 10b). The strengths of the 
vortex filaments adjust to accommodate the changing circulation distribution so that the 
wake always contains the appropriate amount of circulation. EHPIC/HERO also allows 
the option to have filament release points move to the centroids of the zones they 
represent automatically. This usually leads to a physically consistent model once the 
solution converges. Allowing the filaments to move can slow the optimization process 
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and is often not even necessary for this method because the filament release points will 
always be near the centroids of the zones they represent. 

In the second approach, illustrated in Figure 2-1 1, the user assigns fixed values for 
the fraction of the peak circulation to be associated with each zone, and therefore, to be 
assigned to the corresponding filament. Figure 2-11 shows an example where three 
inboard zones (and filaments) are each assigned one third the total circulation trailed 
inboard of the circulation peak and one tip vortex is assigned the total outboard 
circulation. It is usually wise to have the filament release points move automatically to 
the centroid of the zones for this method because the zone boundaries are no longer tied 
to the filament release points. If the zones and Filaments have to traverse long distances 
during the calculation or if they tend to bunch together (as in Figure 2-1 lb), this approach 
can suffer from convergence difficulties. However, when successful, it offers a method 
by which the code creates a physically consistent, evolving wake model that requires no 
post-calculation iterations by the user. Again, the optimization process is delayed when 
filament release points move, but once the calculation is near the optimum circulation 
distribution, the release points will cease to change and the optimization algorithm will 
proceed even more rapidly than the previous method because the fractional strengths of 
the vortex filaments are fixed from step to step. The choice of the best method often 
depends on the user's preference and the particular problem being studied. Additional 
discussion on this point is contained in Reference 12. 


2.7 Investigation of Close Blade/Vortex Interaction Modeling 

Hover performance is in general highly sensitive to the calculation of aerodynamic 
loading in the blade tip region. The tip loading is often influenced by the strong 
interaction with the vortex from the preceding blade. One topic of research during the 
present effort involved methods for more accurate and efficient treatment of this effect. 

A candidate method for accurate resolution of blade-vortex interaction involves the 
application Analytical/Numerical Matching (ANM) (Refs. 28 and 29). This approach 
includes the development of a new way to handle near-field interactions between the tip 
vortex and the rotor blade tip of the following blade. The problem arises because the 
inflow velocity at the blade and the resulting blade loading, which are rapidly varying in 
space, must be defined accurately to assure coirect performance predictions. To calculate 
the velocity induced at enough points on the blade to assure adequate resolution is 
computationally expensive in the context of a vortex lattice model, and would in general 
require a special high-resolution region of the lifting surface. Furthermore, this region 
would have to be adaptive since the position of the vortex is not known ahead of time. 

ANM involves an approach similar to the method of matched asymptotic 
expansions. This approach uses a low resolution numerical calculation in conjunction 
with an analytical near-field correction. As a result, the lifting surface calculation is 
needed at very few points on the rotor blade, with the high-resolution near field being 
constructed from the local analytical solution. To implement this approach, the 
numerical free wake velocity field must be smoothed with an artificially fat vortex core 
when velocities on the rotor blade are computed. Because this smoothing produces very 
gradual variations in velocity, even due to near-field contributions, only relatively few 
control points on the blade are required to reconstruct this velocity field accurately. Note 
that the far-field which lies outside the fat vortex core, and which actually includes most 
of the free wake, is relatively smooth anyway. The fat core smoothing is used only to 
calculate wake effects on the rotor blade, whereas the actual core is used when 
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calculating velocities on the wake itself, so that the vortex filament positions are still 
being accurately computed. 

To implement this method it is necessary to develop an accurate near-field solution, 
ideally an analytical solution based on the local filament configuration. This solution 
should incorporate the position and orientation of the vortex filament passing near an 
idealized semi-infinite blade. Actually, two such filaments must be superimposed. One 
filament adds the contribution of a vortex with a physically realistic core, and the other 
subtracts a vortex filament with the same fat core used in the numerical calculation. The 
net effect in the near field is to cancel the numerical fat core effect and add the effect of 
the actual core size. At the same time, the far-field effect remains unchanged since the 
two portions of the analytical solution cancel in the far field. The vortex filament 
trajectories arc obtained from the numerical free wake calculation modified by the 
analytical solution effect on the vortex. 

Figure 2-12 shows how the blade- vortex interaction can be decomposed into a low 
resolution numerical solution and a high-resolution analytical solution. The numerical 
solution encompasses the full complexity of the problem, except that the fat core 
smoothes out the strong gradients due to the local blade-vortex interaction. The portions 
of the numerical solution outside the fat core radius give the correct velocity contribution 
on the blade. Because of the local smoothing, a relatively low density of vortex 
quadrilaterals can be used on the rotor blade. Basically, the fat core size must be 
comparable to the quadrilateral spacing to assure accuracy. Taken alone, however, this 
approach will give an accurate answer to the wrong problem. The local analytical 
solution then corrects the numerical solution to obtain the correct solution to the original 
problem. 

Figure 2-12 also illustrates the role of the local numerical solution. This solution 
can be obtained from the superposition of a vortex with the actual core size and an 
opposite sign fat core vortex. These vortices can be modeled as straight vortices of 
infinite extent, oblique to the blade. The vortices should be positioned to lie tangent to 
the actual curved tip vortex filament at the point of closest passage beneath the blade. It 
does not matter that these two vortices are straight since they cancel each other in their far 
field, namely at distances beyond the fat core size. However, local curvature effects can 
be added as a refinement, if this is found to be necessary; distant curvature effects are 
handled by the numerical solution. In the near-field region, the opposite sign fat core 
analytical solution cancels the fat core numerical solution, leaving the actual core 
analytical solution. By appropriate choice of the fat core size relative to other problem 
length scales a mathematical overlap region is created, producing a uniformly valid result 
when all the parts of the problem are added together. 

In the course of this effort, the groundwork for the ANM implementation was laid 
down by incorporation of a feature permitting the use of the artificial fat core for wake- 
on-blade interactions. The critical feature needed for full implementation, of course, is 
the near-field solution that characterizes the close interaction of the vortex filament with 
the blade surface (except for the case of direct impingement, which is beyond the scope 
of this study). This near-field solution is composed of two parts; a refined model of the 
vortex core structure and an efficient, high-resolution model of the lifting blade surface. 

A vortex core model was identified for use in this context, drawn from the work of 
Bliss (Ref. 30). The roll-up calculation and core structure model identified in Reference 
30 provides for a three-layer representation of the vortex core, involving a laminar sub- 
core, a potentially turbulent intermediate region, and an in viscid outer roll-up region. 
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Actual Blade/Vortex Interaction 


EQUALS 



Low Resolution Fat Core Numerical Solution 



High Resolution Local Analytical Solution 


Figure 2-12. Decomposition of the actual blade-vortex interaction problem into a 
° low resolution numerical solution and a high resolution local 

analytical solution. 



The analysis in Reference 30 also provides a way to relate some of the important 
integrated properties of the core to spanwise lift and drag loading on the generating blade. 

Identification of an appropriate analytical lifting surface model proved more 
difficult. The approaches presently in the literature (e.g., Ref. 31) involve infinite-span 
lifting wings, while proper treatment of tip effects requires a semi-infinite span to be 
used. Reviews of the available literature on analytical lifting surface methods (e.g., Ref. 
32) did not yield a method sufficiently flexible to be easily incorporated in a local 
solution. Further investigation indicated that a numerical inner solution was called for. 
Previous experience in adapting relatively simple numerical methods to serve as the inner 
solution in an ANM (or 'NNM') implementation indicated that while not ideal, high 
accuracy could nonetheless be achieved while still retaining a net gain in computational 
efficiency. 

Two numerical methods were studied. The first involved the use of a high- 
resolution panel solution for the lifting blade in the vicinity of the tip, using an extension 
of the fixed wing compressible panel method described by Magnus, et. al. (Ref. 33). This 
method, which is implemented in the commercial panel code PANAIR, employs a mix of 
source and doublet singularities to capture both the thickness and lifting effects. 
Calculations of subsonic flow around thick, lifting wings and pressure correlations with 
existing airfoil section data were encouraging. However, it was found that the 
computational demands of the panel method inner solution, even when restricted to an 
isolated high-resolution region close to the tip, has so far kept this implementation from 
being cost effective. 

Investigation into another more promising possibility began recently. Reference 34 
describes an exceptionally simple (hence efficient) numerical lifting surface method 
originally developed for the analysis of unsteady flow over wings with control surface 
deflection. This approach involves distributing a set of point doublet singularities over a 
set of panels on the mean camber line of the lifting surface. Though this method lacks a 
representation of thickness, when reduced to the steady case, the point doublet method 
offers a very efficient method for resolving the loading near blade tips and during close 
vortex interactions. Also, the simplicity and efficiency of the model appears to make it a 
suitable candidate for implementation in the ANM framework presently in place. Work 
on this topic is ongoing, and this refined tip model is a candidate for implementation in 
future versions of EHPIC/HERO. 
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3. 


HIGH RESOLUTION ROTOR WAKE COMPUTATIONS 


The work described in this section was motivated by the high sensitivity of rotor 
performance to the wake structure trailing from the blade tip. As discussed in Reference 
9, EHPIC calculations typically involve using a single vortex to model the wake of the tip 
region. Both experiment and experience with numerical performance correlations to date 
support the use of this approach in many physically important cases. However, it is not a 
fully general representation. This section describes the formulation of a high-resolution 
numerical model of the blade wake designed to be more generally applicable and, in 
particular, to offer greater flexibility in modeling the wake trailed from the blade tip 
region. 

3.1 Rotor Tip Vortex Roll-Up 

Recently there has been much interest in the rotorcraft community concerning the 
flow near the tip of a helicopter rotor blade; much of this has been inspired by the 
possibilities of improved aerodynamic performance resulting from advanced blade tip 
design, e.g., swept-tapered planforms such as the S-76, the UH-60, and the BERP tip. 
The advent of tiltrotor configurations with highly twisted blades that depart significantly 
from conventional design also calls for better resolution of the dominant physical 
phenomena which are present in the rotor wake. Under typical operating conditions, 
these rotor blades all tend to exhibit bound circulation distributions that depart from the 
pattern characteristic to conventional untapered planforms with moderate twist, i.e., an 
abrupt drop in bound circulation near the tip (Fig. 3-1). The presence of gradual tip load 
roll-off (e.g., Fig. 3-2) may lead to incomplete tip vortex roll-up by the time of the first 
blade encounter. This introduces a requirement for an analysis that explicitly computes 
tip vortex roll-up. 

General methods for the analysis of tip vortex formation in rotorcraft are not 
available, though experimental work has yielded significant insights into the structure of 
tip vortices. References 35 and 36, among others, have carried out measurements of swirl 
and axial velocities in rotor vortex cores with the objective of characterizing the core 
structure. This work has also provided empirical support for models of the flow field 
inside the vortex core used in previous calculations . Indirect evidence of the behavior of 
rotor blade wakes can also be obtained from smoke and shadowgraph visualizations 
earned out over the last twenty years (Refs. 1, 2 and 37). These studies have confirmed 
the basic wake structure originally observed by Gray in 1956 (Ref. 38), i.e., a single 
strong tip vortex accompanied by a more diffuse inboard wake sheet. The formation of 
this tip vortex is caused by the rapid roll-up of the wake immediately downstream of the 
rotor tip, which in turn is driven by the large gradient of bound circulation near the tip. 
The presence of such loading distributions on conventional planforms is a well- 
established fact, but significant exceptions to this pattern exist on tiltrotor blades as noted 
above. At present, there exists very little measured spanwise load distribution data on tilt 
rotors (Ref. 39 is one of the few examples), but what docs exist - coupled with available 
computations like that shown in Figure 3-2 - indicates that load distributions have very 
broad peaks located well inboard of the tip and characterized by gradual load roll-off 
outboard. 

Such bound circulation distributions will inevitably lead to more gradual roll-up 
of the trailing vorticity, with the likelihood that not all of the maxim um bound circulation 
will be found in a concentrated tip vortex. Shadowgraph visualizations of model tiltrotor 
wakes have indicated the presence of distinct tip vortices (Ref. 37), though in the absence 
of simultaneous measurements of the circulation of these vortices, it is difficult to reach a 
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RADIAL POSITION (FT) 

Figure 3-1. Typical bound circulation distribution on a representative low-twist 
rotor in hover (CH-47B at thrust coefficient .0068 ). 
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Figure 3-2. Computed bound circulation on a high-thrust tilt rotor blade 
(XV-15/ATB of thrust coefficient .0118 ). 





conclusion as to the fraction of the peak bound circulation they represent. Scaling 
calculations that borrow from work on fixed wing wake roll-up (e.g., Ref. 40) can be used 
to estimate the degree of concentration of the trailed vorticity at the time of first blade 
encounter. Figure 3-3 shows a schematic of the bound circulation on a fixed wing along 
with the trailing wake distribution, as well as the major governing equations of the 
Betz/Donaldson vortex roll-up model described in Reference 40. The model yields an 
estimate of the core radius and circulation strength of the tip vortex of a wing as a 
function of the downstream distance, given the peak value of the bound circulation as 
well as its functional distribution. 

Figure 3-4 shows the evolution of the core properties using this model as a 
function of nondimensional downstream distance for elliptic, parabolic, and linear bound 
circulation distributions. The downstream distance on the horizontal axis is 
nondimensionalized as a function of aspect ratio AR, the wing lift coefficient Cl, and 
distance from the location of the peak circulation out to the centroid of the trailing wake. 
This analysis is not directly applicable to rotary wing tip vortex roll-up because of the 
absence of the symmetry of the trailing wake about the peak circulation as well as 
because of the effects of rotation. It was nonetheless judged to provide a useful guide for 
scaling tip vortex core sizes and for estimating the degree of roll-up that occurs in typical 
rotary wing calculations. 

To investigate this point, assume that the To specified in these figures is treated as 
the peak bound circulation near the tip and that the rotor blade tip wake rolls up into a 
single vortex. If the further suppositions are made that the distance from the position of 
the peak circulation to the blade tip can be equated to the wing semispan in the analysis 
of Reference 40 and that the downstream distance is roughly equivalent to the curved 
path described by the circular wake, then an estimate of the scale of the vortex core size 
and circulation as a function of azimuth angle downstream can be computed. Denoting 
the nondimensional distance from the peak circulation to the blade tip as kgR and 

nor mali z in g the peak circulation by QR 2 yields a measure of the distance downstream of 
the blade for use in the plot shown in Figure 3-4: 

(3-1) 

2 k R 2 V 

(Here, \j r is the azimuth angle downstream of the generating blade in radians). This 
formula applies for both linear and parabolic distributions and can be used as a rough 
guide to the point at which the horizontal axis of Figure 3-4 should be entered to estimate 
die fraction of roll-up that has taken place. For example, on an XV-15 rotor, a typical 

value for r m „ is .03, and the location of most interest is the first blade encounter which 

occurs at v = 2x/3. Assuming that the tip roll-off is parabolic and that kR = .05 (i.e., the 
peak circulation is five percent inboard of the tip), D is 12.56 at the first blade encounter, 
the plot of Figure 3-4 indicates that essentially all of the trailing vorticity will have rolled 

up into the vortex by this point, since the fraction T/To plateaus at 1.0 at values of D 
between 3 and 6. However, for typical XV-15 cases kR is at least 0.15, in which case D 
would be reduced to 1.4. As indicated in Figure 3-4, fora parabolic roll-off the roll-up 
fraction at first blade encounter would then be about 0.8 for parabolic tip loading and 
about 0.5 for linear loading. Realistic loading is probably closer to the former rather than 
the latter, but in either case it seems likely that the tip vortex structure that the following 
blade encounters will not be dominated by a single strong, 
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Aircraft centerline 



Figure 3-3. Schematic of the roll-up of bound circulation from wing with bound 
circulation distribution T(y) on the semispan (Ref. 40). 



Figure 3-4. Estimated trailing vortex circulation strength and radius as a function 
of downstream distance (Ref. 40). 
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fully rolled-up tip vortex. Most of the vorticity will have rolled up, but a substantial 
fraction will remain outside the core, leaving a partially rolled-up filament whose strength 
is substantially below the full value of the bound circulation peak. 

3.2 Structure of Vortex Roll-Up Calculations 

In principle, a variety of options exist for obtaining improved computational 
resolution of the rolling-up vortex wakes. Recent work in this area has included attempts 
to apply both Euler and Navier-Stokes solvers to the computation of wingtip roll-up 
(Refs. 41 and 42). To date, the interrelated problems of high numerical diffusion and 
large computation time have precluded practical application of such tools to the rotary 
wing problem. In the context of EHPIC/HERO, it was judged most appropriate to 
develop a refinement of the present Lagrangian free wake model using vortex filaments. 

3.2.1 Initial resolution of rolling-up near wake 

. In the EHPIC code, an initial guess of the vortex wake geometry is successively 
refined by an iterative relaxation approach based on the influence coefficient method 
described in Section 2. Formally, this iterative relaxation can be considered a systematic 
method for reducing the error in the position of the initial guess relative to the converged 
solution. Within the iteration, various 'error modes' associated with the initial solution 
are successively damped and the solution is marched toward the converged state. 
Typically, several relatively widely spaced filaments are used to represent the inboard 
vortex sheet and a single vortex filament is used to represent the strong rolled-up tip 
vortex. Early attempts to use a high -resolution treatment of the vortical flow near the tip 
region in EHPIC encountered numerical difficulties (Ref. 9). This was a direct result of 
the strong amplifications of the high frequency components of the error modes associated 
with the initial guess (which may be far from the converged wake geometry) and the 
highly nonlinear variations in the wake-induced velocity field. 

In the present effort, a new computational technique was developed that builds on 
the existing wake model to permit a high-resolution treatment of the flow downstream of 
the blade tip. The new method is based on the use of a sequence of coarse-to-fine 
Lagrangian grids and is developed in the spirit of a multi-grid method which is widely 
used in finite-difference and finite-element calculations (Reference 43). The key idea 
behind the present technique is the observation that at the outset of the calculation where 
the approximate discrete solution deviates substantially from the actual 'continuous' 
solution (i.e., the converged wake geometry), the error norm of the guessed initial 
solution is large. If this solution is projected onto a fine grid using small elements, the 
high frequency error mode would amplify rapidly because of the non-linearity of the 
problem (the iterative influence coefficient method is based on a linear Newton's method 
and neglects the non-linear terms in the problem). However, if a coarse grid (i.e., a small 
number of vortex filaments each discretized using a few BCVE elements with large arc- 
length) is used, the fast growing high frequency error components are not present Thus, 
the error norm can be represented as a Fourier series comprising of various Fourier modes 
but with the highest frequency limited by that resolvable on the given grid. 

The slow growing low frequency modes are damped in the iterative relaxation 
process and the calculation converges to an approximate coarse grid discrete solution, 
which is presumably 'closer' to the actual continuous solution than the initial guess. Thus, 
the first approximation to the desired high-resolution solution is thus the single-tip- 
filament solution obtained in previous EHPIC calculations. The error norm associated 
with the coarse discrete single-filament solution is correspondingly lower than the initial 
guess. This solution can be projected onto a fine grid, and because the amplitude of the 
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resulting high frequency error mode is much lower, the non-linearity of the problem is 
not significant (any non-linear term is small) and the discrete solution can be marched 
toward a converged state on the finer grid. The fine grid discrete solution will represent a 
better approximation of the actual continuous solution compared to the coarse grid 
representation. 

Based on the this general numerical strategy, we have incorporated an optional 
algorithm for the high-resolution calculation of flow near a blade tip into the 
EHPIC/HERO code. In this algorithm, a baseline calculation using a coarse-resolution 
wake representation comprised of a single tip vortex filament and several inboard wake 
filaments is carried out. A representative configuration of the wake is given in Figure 3- 
5. This configuration is first iterated toward a converged state with a single strong tip 
vortex as in the EHPIC code. After the initial convergence, multiple copies (depending 
on the desired level of resolution at the tip) of the spatial geometry of the tip filament are 
made and attached to the trailing edge of die blade at stations close to but inboard of the 
tip. An example of this is shown in Figure 3-6, where a converged tip filament has been 
mapped inboard, resulting in a two-filament representation of the trailing wake at the tip. 
Once this initial configuration is set up, the relaxation solution is recommenced and this 
solution is again iterated toward a high-resolution converged state. Usually, convergence 
is achieved within several steps because the error norm of the guessed solution in the 
fine-grid calculation is much lower than would have been the case had two such closely 
placed filaments been started "from scratch". 

The two-filament case shown in Figure 3-6 is only one example of the possible 
application of this technique. To date, calculations with as many as five trailers from the 
tip vortex roll-up zone have been used. Additional refinement is possible within this 
scheme, however computational constraints still limit full resolution of the roll-up 
process, as will be discussed at the end of this section. 

3.2.2 Amalgamation of multiple tip vortex filaments 

In the context of the present scheme, the vortex sheet trailing from the tip, 
discretized into a relatively diffuse set of vortex filaments, typically rolls up into a tighter 
bundle with the filaments effectively collapsing into a single vortex tube far downstream 
of the generating blade. As discussed in Section 3.1, the speed at which this proceeds 
(i.e., the distance downstream required for the roll-up to be completed) scales with the 
size of the tip vortex zone and the steepness of the bound circulation gradient. In general, 
complete resolution of this process is not feasible, nor is it necessary in most practical 
calculations of interest in the current context. A significant question then arises as to how 
to compute the roll-up in a way that leads to robust convergence properties without 
compromising reasonable physical accuracy. In terms of quantities of interest, such as 
surface pressure distribution on the blade immediately following the generating blade, the 
resolution of the rolled-up filaments once they pass by is inconsequential and the multiple 
filaments can be replaced by a single filament while conserving circulation. This is also 
desirable because it reduces the computation time without significantly compromising 
numerical accuracy. 

In EHPIC/HERO, then, the second stage of the high-resolution tip calculation 
involves a merging algorithm that amalgamates the vortex trailers and replaces the 
portion of the multiple tip filaments which are beyond the first blade encounter by a 
single filament. In this part of the calculation, a cross-flow analogy is used to carry out 
the final stages of the amalgamation. The general three-dimensional roll-up problem is 
steady in the rotating reference frame, and thus can be approximated as a two- 
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Figure 3-5. Representative wake configuration of a baseline calculation with a 
single tip filament. 



Figure 3-6. Representative wake configuration of a fine grid calculation starting 
from a converged baseline calculation; tip filament is mapped 
inboard to give a multiple filament representation of the rolled-up 
wake. 



Figure 3-7. Schematic illustrating the merging of multiple tip filaments 
computed using a two-dimensional crossflow analogy. 


28 



dimensional unsteady problem, with distance along the curvilinear trajectory of the 
centroid of the vortex bundle being the 'time-like' variable. The use of a cross-flow 
analogy in three-dimensional steady flow calculations is not new and has been 
extensively exploited in the calculations of separated flow over simple bodies of 
revolution (e.g., missiles) at high angle of attack (Ref. 44). 

This amalgamation procedure requires that the user select an azimuthal station for 
the beginning of the merging process. This station should be well downstream of both 
the generating blade tip and the first blade encounter, a location halfway between the first 
and second blade encounters is typically appropriate. Its exact location is computed as 
the centroid of the multiple free filaments in the initial roll-up at this azimuthal station. A 
trajectory starting from this point and consisting of stations along each of the multiple 
filaments is then computed, with the centroid of the bundle being computed at each 
azimuthal station. A cross-flow plane is then created at each station with its unit normal 
vector computed from the mean tangential vectors of the multiple tip filaments at the 
given station (see Figure 3-7). The multiple tip filaments are then projected onto the 
plane and represented as two-dimensional vortices. 

The time evolution of these quasi-2D vortices is computed in the cross-flow plane 
by integrating the two-dimensional vorticity transport equation. A small time step size is 
chosen such that several steps are executed between a given pair of stations, and a fourth- 
order Runge-Kutta integration scheme is used to ensure high numerical accuracy. At 
each station, the computed locations of the two-dimensional vortices give the cross-plane 
location of the nodes of the vortex filaments. A separation criterion is specified in terms 
of a fraction of the core size of the vortices. If the separation distance between any pair 
of two-dimensional vortices at a given station is less than the separation criterion, the two 
vortices are lumped into a single vortex located at their mutual centroid. The 
corresponding pair of vortex filaments are also merged. Typically, merging of all the 
vortices is completed within five to ten azimuthal stations. 

3.3 Sample Calculation 

The high-resolution tip vortex calculation, coupled with the two-dimensional 
analogy merging algorithm, have been successfully implemented in EHPIC/HERO and 
several test calculations have been carried out The focus of interest here is on rotors 
whose bound circulation distributions exhibit relatively slow roll-up of the tip vortex. 
Test computations have been carried out on an S-76A main rotor as well as on various 
tiltrotor configurations. Selected tiltrotor results are presented here since they best 
illustrate the high-resolution model at work and are in fact potentially one of the most 
important applications for this capability. 

As stated above, computed results on tiltrotor blades often show low bound 
circulation gradients near the tip, indicating that the tip vortex may be only partially 
rolled up by the time it reaches the following blade. Though a concentrated tip vortex 
almost always is present at the first blade encounter, the strength of this vortex tube may 
well be substantially less than the peak bound circulation, and the distribution of voiticity 
in the tube may be poorly represented by a single vortex. To investigate this possibility, 
successively more refined multiple vortex runs were undertaken for the case of an XV- 15 
main rotor. Test data on the hover performance of this configuration is available in 
Reference 45. The blade is 3.81 m (12.5 ft) in radius, with a constant chord of 0.354m 
(1.16 ft.). The tip speed is 234 m/s (769 fps), and the blade has roughly 37 degrees of 
washout, distributed in a nonlinear fashion across the span. Detailed performance 
correlation for this rotor is discussed in Section 6. 
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A baseline calculation was run with the root pitch at 42.5 deg., initially employing 
a very coarse wake with four free filaments trailing from the span: three distributed 
evenly inboard, with a single vortex trailing from the tip (Fig. 3-8). This coarse initial 
condition was chosen to dramatize the effect of the refinement afforded by the multiple 
filament roll-up calculation. This very simplified wake model yields a predicted thrust 
coefficient of 0.0120 and a Figure of Merit of 0.741, substantially lower than the value of 
roughly 0.78 that is estimated from the measured performance data in Reference 45. 
Note that to facilitate this demonstration, the tip filament was trailed from the tip of the 
blade, one quadrilateral outboard of the centroid of the tip circulation zone. This partially 
accounts for the under prediction of the Figure of Merit 

To improve the resolution of the tip wake, three additional filaments were added 
to the high-resolution zone, which extended 75 deg. beyond the first blade encounter. 
Top and oblique views of the geometry of the refined wake are shown in Figures 3-9 and 
3-10, indicating the trajectories of the filaments in the high-resolution region. Each 
filament uses 14 free elements subtending 10 degrees of arc in this region. An estimate 
based on a computation of D in Equation 3-1 suggests that the roll-up should be complete 
in roughly 150 degrees of arc, so elements of this arc length are sufficient to resolve at 
least the full rotation of the tip vortex bundle. The prescribed amalgamation zone applies 
the quasi-2D flow described in the previous section over 70 degrees of arc, leading to the 
smooth merging seen in Figure 3-9. 

The bunching of the three outer vortices, representing 87% of the peak circulation 
strength, is evident in the figure, and allows a more realistic representation of the flow 
field in the vicinity of the blade. The gross features of the roll-up process are captured, 
and the improved resolution is reflected in the improvement in the predicted performance; 
a Figure of Merit of 0.772 at a thrust coefficient of 0.0121 is predicted, close enough to 
the measured performance curve to be within the scatter of the data. Figure 3-11 shows a 
comparison of the predicted thrust distribution along the span, indicating the sharper peak 
that is associated with the implementation of the high-resolution wake model. 

As indicated by these calculations, the high-resolution multiple vortex model 
presently in place has a substantial capability to obtain refined computations of the wake 
trailing from the tip region. However, certain computational limitations still restrict the 
applicability of the present model. The trailing wake should in principle be represented 
by a continuum of trailing filaments rather than by the relatively coarse discretization 
used here. Adding filaments will drive up the CPU time, but it is more the issue of 
robustness that is of most concern. At present, selecting the vortex core radius to be 1.0 
to 2.0 times the distance to adjacent filaments in the high-resolution region produces 
generally consistent, well-behaved results. However, additional spanwise refinement 
brings a requirement to introduce still smaller azimuthal segments to properly resolve the 
wake-on-wake interactions. This leads in general to greater difficulty with convergence, 
though this situation can be aided substantially if the limitations on wake-on-wake 
velocity discussed in Reference 12 (see the input parameter ICON) are invoked. For 
highly refined wakes, though, the overlapping of filament cores becomes inadequate, and 
a higher order model should be used, possibly including vortex sheet elements. 

In addition, this type of inviscid roll-up model is appropriate only to resolve the 
larger scales of the vortex formation process. Refinement of the order of less than a few 
percent of core radius requires implementation of a direct model of viscous and turbulent 
flow within the core, as discussed in Section 2.7. To date, adding this level of refinement 
has not proved cost-effective for the performance optimization applications explored 
here. 
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Figure 3-10. Side/oblique view of the converged high-resolution solution for the 
XV- 15 rotor (wake of one blade only shown). 



Radial distance (ft) 


Figure 3-11. Bound circulation distribution with initial coarse-resolurion and 
Final high-resolution wake models for the XV-15 case. 
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3.4 Side-edge Separated Flow Modeling 

In the present effort, an approximate method for computing the effect of flow 
separation on the side-edge of a blade has been developed. In general, this is a viscous 
phenomenon whose details cannot be computed with the present in viscid model. The 
motivation for this particular feature of the model was to attempt to develop a simple 
model that quantifies the likely effect of the lift-off of tip vortices from the blade edge, 
and to permit the user to assess its importance in particular cases. 

The present method represents an extension of the two-dimensional cross-flow 
calculation developed for the merging of multiple vortex filaments, as described in the 
previous section, to the edge flow problem. In this calculation, the cross-flow analogy is 
applied to the treatment of the flow along the side-edge of the blade, and the three- 
dimensional steady separated flow is computed as a two-dimensional unsteady separated 
flow problem. In this case, cross-flow planes are computed which are in-plane with the 
unit normal vector n of the blade surface and the radial direction of the blade (Fig. 3-12). 
The plane is placed along the edge of the blade and the flow velocity normal to the blade, 
resulting from the angle of attack of the blade airfoil, is projected onto the plane. This 
results in a two-dimensional unsteady problem involving a flow past a flat plate normal to 
the flow direction. Time integration of this two-dimensional problem gives the location 
of the separated vortex along the side edge of the blade. At t= c/Uoo, where c is the chord 
length of the blade and Uoo is the projected free-stream velocity in the two-dimensional 
case, the two-dimensional calculation gives the location of the separated vortex near the 
trailing edge of the blade, which is used as an attachment point for the tip vortex in the 
three-dimensional case. 

Unsteady separated flow past a normal flat plate is a classical problem that has 
been extensively studied (Refs. 46 - 48). The work of Pullin (Ref. 48) is of particular 
interest and can be directly applied to the present calculation. In his model, the self- 
similar rolled-up vortex sheet emanating from an impulsively started flat plate is 
computed. Based on his self-similar results, it is straightforward to derive the following 
relation between the location of the separated vortex and the total circulation in the vortex 
and time: 


Z = 


4 2.64 


1/2 


(3-2) 


where Z is the height of the vortex above the blade, T is the circulation of the vortex and t 
is time. This relation suggests that given a vortex with known location and circulation, 
the locations of all other vortices with different circulations can be estimated based on 
this similarity. 

This has been tested by computing the separated vortex of a two-dimensional 
normal flat plate. Two separate flow configurations were examined, which resulted in 
two vortices of different net circulations.^ It was found, however, when the coordinates 
are properly scaled, the two vortices appear to be similar, with nearly identical locations. 
This is shown in Figure 3- 13a and 3- 13b and the similarity argument is well proven. This 
similarity relation considerably simplifies the implementation of a simple side-edge 
separation model in EHPIC/HERO since we need only to compute the location of the 
vortex once (at the first step for a given blade pitch angle); at all subsequent iterations, 
the location is obtained by using the scaling in Equation 3-2. The effect of this model 
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Figure 3-12. Schematic showing blade and cross flow plane at side edge of blade. 



a) Normalized circulation = 1.0 



b) Normalized circulation = 2.0 


Figure 3-13. Demonstration calculation of a self-similar 2D unsteady vortex 
solution, modeling the flow development on a blade edge. 
Calculations are at equal times. 
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when invoked is to alter the vertical release point of the tip filament. The inputs required 
for this model are described in Reference 12. 
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4. FORMULATION OF THE ROTOR BLADE STRUCTURAL MODEL 


The structural properties of the helicopter blade are in many cases important to the 
evaluation of the performance of a helicopter in forward flight Thus, one of the major 
efforts within the development of the EHPIC/HERO code has been to incorporate a 
realistic finite element (F.E.) representation of the blade. This section discusses the 
formulation of the structural model together with its capabilities and limitations and the 
manner in which it is coupled to the aerodynamic wake model. Further description of the 
inputs required for this portion of the analysis are given in Reference 12. 


4.1 Finite Element Structural Model of the Helicopter Blade 

The particular finite element (F.E.) model used here to represent the helicopter blade 
accounts for extension, twist and transverse bending displacements. To accurately 
simulate these deformations, the blade is discretized into a number of beam finite 
elements each having a total of 14 degrees of freedom (d.o.f.). Stiffness properties for 
each element are computed from the cross-section geometry and material properties 
supplied by the user. Similarly, the blade mass distribution is used to both define the 
nodal forces due to blade rotation and also the contributions of blade rotation inertia 
forces to the stiffness matrices (geometric stiffening). The resulting elemental stiffness 
matrices are then assembled and any constrained d.o.f. eliminated to finally yield the 
global stiffness matrix for the complete blade structure. The approach taken is similar to 
previous implementations of F£. methods for rotorcraft applications, such as Reference 
49. 


The transfer of information between the structural and the aerodynamic models is 
achieved by requiring force equilibrium and by specifying the kinematic relations that 
define the blade geometry in terms of the undeformed blade shape and the vector of nodal 
deflections. In essence, nodal forces, f, are computed from the distributed aerodynamic 
forces and the blade rotation. The steady state deflections are determined by solving the 
linear system of equations, [K]s=f. Finally, the deformation vector is used in conjunction 
with the undeformed blade geometry to furnish the deformed blade geometry. The 
geometry update perturbs the flow-field which in turn alters the nodal forces. Thus an 
iterative process is invoked until convergence is attained. This iterative process is done 
in parallel with the wake relaxation calculations so that the entire process essentially 
converges to the final deformed blade and associated flow-field. The remainder of this 
section explains in greater detail the derivation of the structural model and the manner in 
which it is coupled to the wake model. 


4.1.1 Assumptions 

The assumptions inherent in the blade model and geometry are stated below; 

• The blade displacements are of sufficiently small magnitude that: 

- linear constitutive relations between stress and strain are applicable, 

- the transformation matrix relating the local axes of each element may be 
regarded as constant and equal to the corresponding matrix in the 
undeformed state. 
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- rotations due to deformation are assumed to commute, and, 

- the twist, bending and extension deformations may be linearly 
superimposed. 

• The blade material is assumed isotropic and the stress-strain relation obeys 
Hooke’s law. 

• The elastic axis for each element is defined. The elastic axes of any two 
adjacent elements coincide at their mutual joining section (see Fig. 4-1). In 
other words, the elastic axis is continuous along the blade. This is necessary to 
correctly define the assembly of the individual blade F.E. s. 

• The principal axes of the cross-section for each element are assumed to be 
perpendicular to the elastic axis of that element. This implies that if there are 
sweep and anhedral changes between consecutive elements then their principal 
axes will not coincide at their mutual section. The degree of approximation 
introduced into the bending calculation will increase with the amount of sweep 
and anhedral change between adjacent elements. It will also decrease with 
slenderness of the element since the discrepancies resulting from non-alignment 
of the principal axes occur locally in the neighborhood of the joining section. 

Note that the last two assumptions are mainly due to the fact that warping effects are 
modeled in the analysis. One of the chief advantages of the finite element method is its 
versatility in the assembly of the constituent elements. For simple elements, e.g., pure 
beam elements and bar elements, one is free to assemble the components in whatever 
orientations one chooses. Furthermore, discontinuities in the mass and stiffness 
properties from element to element are permitted. However, when modeling warping 
deformations the line of shear centers, or elastic axis, plays a significant role. The current 
formulation approximates the elastic axis by a sequence of straight line segments and it is 
the desire to accurately represent the elastic axis that results in the preceding last two 
assumptions. Thus to the extent that warping effects are significant, failure to satisfy the 
last two assumptions and suitably approximate the elastic axis leads to error in the 
solution. In most cases however, and for the closed tubes representative of helicopter 
rotor blades, warping effects will be dominated by deformations arising from pure 
bending and torsion, and thus violation of these assumptions will not lead to significant 
error. This has been verified by numerical testing of the F.E. model for loaded structures 
containing 90° elbow joints and discontinuities in the beam stiffness properties. 


4.1.2 Blade geometry 

Each of the blade segments defined in the EHPIC/HERO blade geometry input 
corresponds to a single structural finite element This is assumed to be adequate for the 
hover situation where static deflections are sought The global axes for the assembled 
blade are denoted by XYZ corresponding to the blade axes defined in Figure 2-5. Local 
axes, xyz, are defined for each element such that the x-axis coincides with the elastic axis, 
or the line of shear centers, of the element. Axes y and z are derived from the 
transformation applied to the global Y and Z axes as described below. The 
transformation matrix relating the local element axes to the global axes is derived from 
the local segment layout specifications. The segment geometry is specified as follows 
(see Fig. 4-2): 
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Figure 4-1. Correct and incorrect alignment of the elastic axis (E.A.) 
between adjacent elements. 
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(1) The planform is first defined. Each blade segment has length, SL, along the 
global X-direction and chord length, c, in the global Y-direction. The sweep, A, defines 
orientation of the quarter chord line for the segment Note that for non-zero sweep, the 
length of the finite element along the quarter chord length differs from the length 
measured along the blade X-axis. If one finite element is associated with each blade 
segment then the element length shall in fact be: 


1 = 


SL 

cosy 


(4-1) 


where y is the anhedral (see step 3 below). 

(2) A camber and then a pre-deformation twist gradient are defined over each 
segment. This information is not included in the transformation matrix since any effects 
due to camber and pre-twist upon structural properties can be more accurately specified 
in the information on blade cross-section properties (see Reference 12). Addition of 
camber would be reflected in the cross-sectional moments of area and pre-twisting would 
affect primarily the orientation of the principal axes. These parameters are directly 
specified in the blade cross-section input file discussed in Reference 12. 

(3) Anhedral is then applied to each segment about an axis parallel to the global 
Y axis and passing through the left hand end (nearest to the rotor hub) of the segment. 

The direction of this rotation is in the negative Y-direction., i.e., positive anhedral, y, 
results in the blade drooping down. 

(4) Finally, collective pitch in the form of a rotation about the global X-axis is 
applied to the assembled structure. 

This sequence of rotations is used to define the transformation matrix relating the 
local axes to the global ones of the EHPIC/HERO code. An additional 180° rotation 
about the global X-axis precedes the above rotations since the local finite element z-axis 
is positive upward whereas the global Z-axis is positive in the downward direction. The 
preceding parameters are supplied in the blade geometry input file. 

From the above sequence of rotations the local and global axes are related by: 

+ T c T y T A t 180° jyj (4-2) 

where Tjgo 0 , T A , Ty and T c are the transformation matrices corresponding to the 180° 

rotation, sweep, anhedral and collective operations respectively, and the coordinates, 
XqYqZo , are the global coordinates of the origin of the local axes. Here the origin lies 
on the elastic axis at the end of the clement nearest the rotor hub. The combined matrix, 



[ T rot] = T c T y T A T 180° 


c^A c^A Sy 

c c s A' s c s y c A “ c c c A”^c s 'y®A ^c^y 
L S c S A 4C c SyC A -ScC A +C c SyS A -CgCy _] 


(4-3) 
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where sm and C(.) denote sin(*) and cos(») respectively. Note that the rotations due to 
deformation can also be referred to the global axes using this transformation since the 
deformations are assumed small and the rotations thus commute. 


4.1.3 Element degrees of freedom 

The specification of the shape functions and the fourteen degrees of freedom of each 
element is summarized here. Each element has two end nodes and one node at its 
midpoint, as shown in Figure 4-3. The degrees of freedom correspond to translational 
and rotational deformations at these nodes. The deformation of the element at any point 
is estimated by interpolation of the nodal displacements using the shape functions. Let u, 
v, and w denote the displacements along the local x, y, and z axes respectively and let 

0 denote the twist deformation about the x axis. Then the generalized displacement 
vector is defined as: 


Lag: 


Twist: 
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(4-4) 


where the subscript (»), x denotes the derivative with respect to the local x axis 
coordinate and 1 is the element length. Thus qj and q2 refer to the displacement and 
corresponding slope due to bending in the y-direction at the left hand node. The 
corresponding right hand node deformations are q 3 and q 4 , and so forth for the other 
displacements. Note that the slopes, v, x and w, x , can be regarded as a small positive 
rotation about the local z-axis and a small negative rotation about the local y-axis 
respectively. 

The transverse displacements, v and w, are interpolated using cubic Hermitian 
polynomials as is the common practice in beam finite element formulation. Quadratic 
polynomials are used to interpolate the torsional and axial deformations. This is the 
simplest element interpolation scheme yielding a consistent formulation for coupled 
torsion-bending (Reference 49). Specifically: 
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where, {O 3 } = i 


(4-6) 
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and £=x/l . The preceding relations may be expressed in compact form as: 
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(4-7) 


where (q) is the vector of generalized displacements and [ O ] is a (4x14) matrix 
appropriately constructed from Eqs. (4-5) and (4-6). 

Finally in the formulation of the elemental stiffness and mass matrices it is valuable 
to define principal axes of the cross-section, T| and £ , which are oriented such that: 

J TlC dA = 0 (4-8) 


The angle P is then the angle between the local y and tj axes. 


The global degrees of freedom are obtained by resolving the deformations along the 
global axes using the transformation matrix derived previously. At an end node, all of the 
three translational and three rotational d.o.f. are available (since the slopes of the 
transverse bending displacements correspond to rotations). Thus, the translation between 
the local element d.o.f.s and the global ones is achieved using the transformation matrix, 

^rot • 



f Rx l 
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r y 


-w,x • 

l Rz J 

. v *x . 


(4-9) 


where Rx , Ry > and Rz » are rotations due to deformation about the global XYZ axes 
respectively. At the mid-node the preceding translation is more involved since only two 
d.o.f., the twist and axial deformations, are available in the local axes and additional 
constraints are necessary to uniquely define the twist and stretch in the global directions. 
One approach would be to specify the four remaining local d.o.f at the mid-point by 
interpolating from the end-nodes using the element shape functions, i.e. evaluating 
v(l/ 2 ), w(l/2) and their slopes using Eqs. (4-5)-(4-7). These together with the local twist 
and extension deformations completely define the six local displacements from which the 
global deformations readily follow. However, it was found that this led to numerical 
problems in the resulting transformation matrix, since the complete element is singular 
for certain blade geometries. This might be expected from the observation that 18 global 
displacements (6 at each node) have been defined in terms of only 14 element d.o.f. 
Hence, the inverse transformation from global to local deformations is in fact non-unique. 
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The approach taken here is to simply define the global deformations to coincide with the 
respective local ones at the mid- node, i.e., 



(4-10) 


This both simplifies the transformation and results in an orthogonal element 
transformation matrix, i.e., if the elemental transformation which will be composed of 

elements of [r rot ] is denoted by Ugl] so that qg = [TglI Al t^ en ITgiJ^ = UglI^ • 
An alternative procedure would be to eliminate the mid-node d.o.f. using static 
condensation. However, this is unnecessary in light of the small number of d.o.f. s of the 
fully assembled model, the additional programming complexity and the further 
approximation that would thus be introduced. 


4.1.4 Derivation of the element strains and stresses 

In order to compute the elemental stiffness matrix the strains arising from the 
preceding displacements must be evaluated. The nonlinear expressions for the strains are 
stated: 

- u.x + ( '•'9.x ),x +5 ( n 2 + C 2 ) 8.x 2 

1 o 

+ 2 v, x z * v »xx { T| cos(P+0) - £ sin (p+0) ) 

+ \ w, x 2 * w »xx ( C c °s(P40) + H sin(P+0) } (4-1 la) 


Yxri « 2 e XTl = (¥*,-£ )(e* + 0ni) (4-llb) 

YxC = 2e xC = (^< + T l)( e ’X + e nl) (4-llc) 

and all other strains are assumed zero. Here, 0 n j is a nonlinear second order torsion 

term, and v F(x,T|,Q is the Saint Venant warping function expressing the out-of-plane 
displacement, u wai p , due to torsion: 

u waip( x * 7 l«C) = ^(x,Tl ,Q 0»x (4-12) 

The linear expressions are easily obtained from above: 

e xx = u *x + ( TOx )>x 

- v lXX { tj cosP - £ sinP ) - w tXX { £ cosP + 11 sinP } (4-13a) 

Yxri = )e,x ; YcC = (^<+n)e,x (4-i3b,c) 
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These strains are expressed in terms of the vector of generalized d.o.f., {q} , and the 
shape functions and their derivatives wj.t x , by substituting for the occurrences of u, v, 

w, and 0 using the expressions, Eqs. (4-ll)-(4-13). This results in: 


-xx 


f - [ TjcosP - CsinP ] {^>3 } 
_ ^ -[Ccosp+TisinP] {<D 3 } 
'P.x (^2 ) + ^(^2 } 

i®2) 

= { Bi } T {q} 


(q) 





0 'j 

. Tfcnl 

L {'•v? 

L 

0 

> > 

.TxcJ 

I'T.C + t,] 

n 

{^2 ) 

L 0 J 


(q) 


{ {B2) t 

1 {B 3 ) T 


{q} 


(4- 14a) 


(4-14b,c) 


The corresponding stresses are derived from the Hooke’s Law: 


a xx = Ee xx 


(4- 15a) 


°xtl = ; <*xC = Gy x £ 


(4-15b,c) 


4.2 Derivation of the Equations of Static Equilibrium 

The discussion above defines the relationships for stresses and strains for the current 
F.E. formulation. The required equilibrium equations are obtained via the principle of 
virtual displacements where virtual work expressions are constructed by considering the 
internal and applied forces subject to virtual displacements. In this framework, the usual 
material stiffness properties are represented by an internal virtual work expression, W‘, 
which in this case is equivalent to the variation of the strain energy expression. Since the 
internal virtual work can be written down directly, it is unnecessary to execute the 
intermediate step of obtaining an expression for the strain energy. The rotating XYZ 
reference frame gives rise to both inertial forces and geometric stiffening which can both 
be conveniently represented by an external virtual work term, W c . The aerodynamic 
forces computed in the wake analysis constitute a further external virtual work term. 
Finally, the equations for static equilibrium are obtained by equating W‘=W e . 

The structural stiffness matrix K is composed of three terms. Kg, K R , and Kq, 
where Kg is the usual material stiffness matrix obtained from the internal virtual work, 
and Kq and Kr are contributions due to the blade rotation. Kr is due to the centrifugal 
loading. This distributed inertial load gives rise to both a nodal force vector due to blade 
rotation and also a softening of the blade. The latter is physically due to the fact that 
displacement of a blade element places it at a different radial location and therefore alters 
its loading. A simple example of this is the case of deformation along the X-axis, u, at a 
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radial location, r, for a blade rotating with angular velocity, £2. A deformation u increases 
the centripetal acceleration from £2r 2 to £2(r+u) 2 with a proportional increase in loading. 

This load increase induces a further displacement, 8u, and in the linear F.E. model is 
equivalently accounted for by Kr. In the steady state deflection case for hover, Kr is 
obtained by developing an external virtual work term for the virtual displacements in the 
presence of the distributed centrifugal forces. The contribution Kq accounts for the 
geometric stiffening due to the axial tensile loading induced by blade rotation. The 
geometric stiffening plays a major role in the range of angular velocities typical for 
helicopters. This effect is analogous to the buckling problem where the axial force is 
compressive and generates internal distributed moments when the blade deforms. 

For simple blade bending examples the blade stiffening effect, Kq, is found to be 
twice the softening contribution, Kr, leading to a net increase in blade stiffness, as is 
clearly the expected behavior. Kq is derived by expressing the forces and moments due 
to rotation explicitly and regarding these as externally applied loads that are accounted 
for by an additional external virtual work term. The external virtual work expression then 
results from these applied loads undergoing virtual displacements. This approach is 
frequently adopted in buckling analysis and bending and torsion problems where axial 
forces are present 


4.2.1 Construction of Kr 

The internal virtual work due is given by: 

W* = <*xx*8Exx + °xT| * ftficri + • &Y X £ dV (4-16) 

Substituting for the stresses and breaking up the volume integral: 
i 1 r 

w = J ( J Eexx • 8exx + G'fcq • 5y xT | + Gfc$ • dA } dx 

(4-17) 

Substituting for the strains using Eq. (4-14) and performing the integrations results in: 
W i «5a T [K E ) a (4- 18a) 

where the desired stiffness matrix is 


k e = 


! lj E{B 1 }{B 1 } T + G({B 2 }{B 2 } T +{B 3 }{B 3 } T )dA )dx 


(4-18b) 


The construction of Kg requires a sequence of integrations, the first being an area 
integration over the area of cross-section at a given station, x , along the element, and the 
second being the integral along the length of the element Evaluation of the area integral 
results in expressions containing various properties of cross-section multiplied by the 
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shape functions and their derivatives w.r.t x . These properties include the cross- 
sectional area, moments of area, area centroids relative to the q and £ axes, and a total 

of nine integrals involving the warping function, X F. The finite element implementation 
employed in EHPIC/HERO does not compute these properties, but instead requires that 
the various cross-sectional area integrals be input direoly via the blade cross-section 
input file. The finite element code requires that these properties be specified at the end 
nodes of each element and assumes that they vary linearly between the end nodes. The 
final integration along the element length involves products of the shape functions and 
their derivatives and is effected numerically using Gaussian integration. A list of the 
cross-section area integrals required in the analysis is given in Reference 12. 


4.2.2 Construction of Kr and f 01 

In order to determine the contributions to the stiffness matrix and the nodal forces due 
to blade rotation one first defines the position vector of a point on the blade in blade 
coordinates, 

R(X,YZ) = XI+YI+ZR 


or, i 

i><X 

ii 

~hT>< 

o o 

[ + [Trod " 

xie + u 

v+qcos(|3+8)-Csin(P+0) 

► 

(4- 19a) 

l 

ZJ IZqJ 

1 

l w+rj sin(P+0)+£cos(P+0) . 




which in local coordinates is, 
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[ZoJ 

1 

. w+T| sin(P+0)+^cos(P+0) . 


where xj e is the distance along the elastic axis of the element containing the point. Then 
the body force at any point on the blade due to rotation is: 

f = - p (OR) x {(OR) x R } (4-20) 

where p is the density of the blade material, and the unit vectors are aligned with the 
global axes. When expressed in the local coordinate system of a particular element, this 
becomes. 


ffxl 
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r x oi 

■ 1 00‘ 


[ x ie + u 
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[TroJ' 
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l fz J 
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l 0 J 
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l w+T| sin(P+O)+^cos(P+0) . 

J 


(4-21) 


The moment about a point on the element elastic axis, Rc a = Y eal+^eaK > due to 

the rotational force acting on a volume element located at R somewhere on the blade is. 
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ID = (R-Eea) x f-dAdx 


(4-22) 


The net forces and moments at a point on the elastic axis defined by x ea are: 


R 


E = 



I f(x,y,z) dA ) dX 


(4-23a) 



J [R(X,Y,Z) - R ea (X ca ,Y ea ,Z ea )] x f(X,Y,Z)dA } dX 

(4-23b) 


where R is the value of X^, at the blade tip. The domain of integration extends from 
X ea to the blade tip since the net force and moment vectors vanish at the blade tip. 


The external virtual work for the inertial forces is obtained by considering the body 
force due to blade rotation, f , as a distributed external force. Imposing virtual 
displacements upon the blade under this distributed load results in the formation of the 
external work term: 


1 




f x 5u + fy£v + f z 8w + (yf z - zfy)50 dA dxj c 


(4-24) 


The remaining procedure is laborious, but straightforward and is briefly summarized 
below: 


• Resolve Eqs. (4-19) to (4-24), in the local element coordinate system. 

• Substitute for all occunences of u , v , w , and 8 , and their derivatives using 
Eqs. (4-5) and (4-4). 

• Replace sin(|3+0) and cos(P+8) by the small 0 approximations. 

• Substitute for f in Eq. (4-24) and discard all terms of order higher than 2. 

• Carry out the cross-section area integrations. As in the computation for Ke, this 
area integral can be directly expressed in terms of certain cross-section 

properties. Since the blade material density, p, is now present in the analysis 
these cross-section properties will be quantities such as the mass per unit length, 
cross-sectional center of mass, torsional moment of inertia per unit length, etc. 
The complete list of parameters needed is given in Reference 12. 

• Finally, evaluate Eq. (4-24) from x} e =0 to xj c =l using Gaussian integration. 
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The resulting integral assumes the form, 


1 

W£ = - 5q T £2 2 J [A 0 ]a-{B 0 }dx 

or, Wjj = 5 q T ( -[KR] a + fot) 

where the contribution to the stiffness matrix, 
1 

K r = £J 2 f [Aq] dx 


and the nodal force vector due to the centrifugal forces, 
1 

f-tfj {B 0 } dx 


(4-25) 


(4-26) 


(4-27) 


4.2.3 Construction of Kq 

The virtual work expression for the net blade rotation forces undergoing virtual 
deformations that accounts for the geometric stiffening effects is stated (see Reference 
49): 

1 

w£ = J {F x (v'Sv'+wW) + M y S(v"0) + M z S(w H 0) 

+ Qe’8(0') + j M x 5(v"w'-w'V) } dx ic (4-28) 

Here, 

Q = ^i Fx(y2 + z2>dA (4 " 29) 

and the terms, F x , M x , My, and M z are simply the local components of the net forces 
and moments due to rotation summed over the portion of blade lying outboard of the 
point xj e on the elastic axis. Note that the first term in Eq. (4-28) represents the usual 
additional stiffening due to an axial force. The virtual work contribution for this first 
term may be viewed as a differential moment arising from structural deformation, 

dM z = F x (v'dx) and dMy = F x (-w'dx), moving through virtual rotations, 8v' and 
5(-w*). 


Eqs. (4-19)-(4-23) and (4-28)-(4-29) contain all of the information necessary for the 
computation of [ Kg ]. Again the actual computation is laborious, but straightforward. 
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One begins by substituting for f , m x , F x , M, and Q in Eq. (4-28) and discarding all 

terms of order higher than 2. The sin(») and cos(») are replaced by their small 0 angle 

approximations, and the occurrences of u , v , w , and 0 are evaluated from Eqs. (4-4) 
and (4-5). The cross-section area integrations arc performed to define mass moments of 
area. Finally, use of integration by parts where possible simplifies the integration along 
the element from xi e =0 to xj e =l of Eq. (4-28). For example, the first term, 

1 

J" {F x (v'5v'+w'8w') dxj e = 

1 

' x ie 

+ f x { J ( v'Sv'+w'Sw' ) dp} dxj e (4-30) 

0 


x ie 


-i "ie * 


F x (xie) J v'8v'+w'5w'dp 


xi»=0 


The quantity (v'5v'+w'5w') is easily evaluated from Eq. (4-7) in terms of the element 
shape functions and the generalized vector of nodal displacements, q . Thus the integral 
contained in the brackets {•) can be written down analytically. The final integration 
along the element from xj e =0 to xj c =l is done by Gaussian integration, and results in, 

W®=-5fl T [KGla (4-31) 

7 

where the geometric stiffening matrix, Kq, is proportional to Q . 


4.2.4 Computation of I 3 " 0 

The coupling of the structural model to the aerodynamic wake analysis is achieved by 
evaluating the nodal force vector, f 3 ® 1-0 , due to the distributed aerodynamic forces and 
moments. The associated external virtual work term is derived in a very similar 
procedure to that used in developing the rotational force vector, f ot . The wake analysis 
presents the finite element routines with an array, QFRC(iseg,ic,k), which represents the 
aerodynamic force on chordwise panel, ic, of blade segment, iseg, in the global blade axis 
direction, k. The index k=l,2,3 denotes the aerodynamic force in the X, Y, Z directions 
respectively (thus k=2 denotes the drag and k=3 the negative of the thrust for the 
chordwise panel), and k=4 represents the aerodynamic pitching moment about the X-axis. 

To derive the corresponding nodal forces, the aerodynamic forces per unit length is 
first obtained by dividing entries of QFRC by the panel width, dX. The resulting vector 
is rotated into the local reference frame to obtain local forces and moments per unit 
length. 
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(4-32) 
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The external virtual work expression for the aerodynamic forces is then (c.f. Eq. (4-24)) 

1 

= J* f|5u + f^8v + f^8w + m*60 + mj8(-w >x ) + m*6(v, x ) dxj e 
(4-33) 

which is evaluated as before by substituting for the deformations, u, v, w, and 0, and 
numerically integrating to obtain, 

= 5a T f aero (4-34) 


4.2.5 Equations of equilibrium 

The derivation of the virtual work terms is now complete. The equations of motion 
follow immediately from equating the internal and external virtual work terms, W i =W e , 
using Eqs. (4-18), (4-25), (4-31), and (4-34): 

5fl T [K E ] a = 5q T ( -[K R ] a +frot) - 5a T [ K g ] a + Sa 7 ^™ 
or, since the virtual displacements are arbitrary, 

[K E + K R + KG]a = F l + f 3 * 0 (4-35) 


4.3 Assembly of the Global Stiffness Matrix and Forces 

The preceding equation, Eq. (4-35), implies an assembly process of the individual 
element stiffness matrices and nodal forces into a corresponding global stiffness matrix 
and applied force vector for the complete helicopter blade. The assembly process 
involves three sub-procedures: the first involves referring the elemental matrices and 
nodal forces to the global axes, the second defines the array indexing that relates the local 
degrees of freedom for each element to the global ones, and the final step entails 
implementation of the boundary conditions at the blade root 

Rotation of the element matrices and nodal forces into a global coordinate frame is 
accomplished in the standard manner: 

[ Kglobal 3 = [Tgl] t K local 3 PglJ"^ (4-36a) 

(f rot }global= [TgiJ (f 01 ) local (4-36b) 
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(4-36c) 


{ f 86 ™} global = [T G J {f^Jtocal 

as may be easily verified by noting that the potential energy and the virtual work are 
independent of the choice of reference frame. Here, [Tql] is the transformation matrix 
described in Section 4.2 relating local d.o.f., a, and global generalized d.o.f., s: 

£ = [TglJ a- 

The blade elements are then laid end to end in sequence from blade root to blade tip. 
Global deformations are defined as outlined in Eqs. (4-9) and (4-10). However, the 
ordering of the global degrees of freedom is different from the local ones, Eq. (4-4). 
Each element degree of freedom is associated with a global one via an indexing array or 
splay matrix, C(k,ie) , where k is the local degree of freedom (k=l,2,...14) , ie is the 
element number, and C(k,ie) is the global degree of freedom. Having specified a C(k,ie) 
for each element then the global matrices may be constructed by ‘splaying’ components 
of the elemental mass and stiffness matrices into their corresponding positions in the 
global matrices. For example, the [i j] entry of the elemental stiffness matrix for finite 
element, ie, is added to the [ C(i,ie),C(j,ie) ] entry of the global stiffness matrix. In like 
manner, the global nodal force vector is built up from nodal forces for each element. 

It remains to specify the actual ordering of the global degrees of freedom. Degrees of 
freedom are numbered upwards starting at the blade root Using the definitions for global 
displacements and rotations given in Section 4.2, Eqs. (4-9) and (4-10), the global 
degrees of freedom are summarized in Table 4-1 (see also Fig. 4-3). 

The construction of the global stiffness matrix is completed by imposing the 
boundary conditions at the root For articulated blades, it is implicitly assumed in 
EHPIC/HERO that the blade is freely hinged in both flap and lag directions, but that the 
remaining degrees of freedom at the root - the three translational displacements and the 
twist about the X axis - are constrained. The boundary conditions are implemented by 
simply deleting the rows and columns of the global stiffness matrix corresponding to 
these four degrees of freedom. For cantilevered blades, by definition all root 
deformations are zero and thus all six degrees of freedom at the root must be removed. 


4.4 Solution of the Structural Equations 

For a given nodal force vector, it is required to solve the finite element 

equations. 


[K]s = £ (4-37) 


for the global deformation vector, $. The iterative solution process whereby the final 
deformed blade geometry and associated wake structure are obtained, gives rise to a 
series of nodal force vectors, £, which changes as the flow-field is updated. For a fixed 
undeformed blade planform [K] remains unchanged so that Eq. (4-37) is repeatedly 
solved for a sequence of right hand sides. The LU decomposition procedure is most 
amenable to this type of problem since the LU decomposition of [KJ need only be 
performed once. The solution £ for any f is then obtained by a computationally cheap 
backsubstitution procedure. 

The end-to-end layout of the beam elements and numbering of the d.o.f. results in a 
banded stiffness matrix with maximum bandwidth, mb=22. In EHPIC/HERO, the non- 
zero elements of [K] are stored in an array B according to. 
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(4-38) 


B(j-i+14 , i) = [KJy 

This is the storage form required by the pair of efficient banded LU decomposition and 
backsubstitution subroutines DECB and SOLB of Reference 50. The LU decomposition 
is performed once by calling DECB immediately subsequent to the construction of the 
array, B. The nodal deflections are then obtained for each of the nodal force vectors by 
the backsubstitution routine SOLB. 

For certain wake layouts and blades, the iterative process has been observed to 
become unstable, particularly for high resolution wake models or relatively soft blades. 
This is believed to be due to the change in structural deformation, and hence the blade 
geometry, that takes place between consecutive iterations. If this change is too large, 
kinks in the trailing filaments develop and undermines convergence. To counter this 

possibility, the maximum change in deformation state during an iteration. As, is limited 
by the parameters DLD and DLR supplied by the user. Here DLD is the maximum 
change in displacement deformation between EHPIC/HERO iteration steps, and DLR the 

corresponding ma ximum change in rotation deformation. If any element of As is found to 
exceed these limits, then the whole vector, As, is scaled so that the maximum change 
requirement is satisfied. The deformation state is then updated, s^-fi+As, using the scaled 
version of As. This limiting is not applied during the initial relaxation since this portion 
of the calculation tends to be relatively insensitive to characteristic changes As. 

Since the overall process is iterative, a tolerance parameter, analogous to those used 
to test for wake convergence, is supplied to verify that the structural computation has 
converged. The required test for all i takes the form: 

displacement deformations: Asj < CONVGS * R (4-39a) 

rotation deformations: Asj < CONVGS (4-39b) 

where CONVGS is the user supplied tolerance parameter. 


4.5 Update of Blade Geometry 

The coupling of the structural model is completed by specifying the relation between 
the global deformation vector, s, and the geometrical parameters used to specify the blade 
layout As described in Reference 12, the blade geometry is defined by the variables: 


SL(iseg) 

SWEEP(iseg) 

TWG(iseg) 

TWR(iseg) 

ANH(iseg) 

CHORD(iseg) 

COLL 


segment X-length 

segment sweep 

segment twist gradient 

twist at left hand end of segment 

segment anhedral 

chord length at left hand end of segment 
blade collective. 


In EHPIC/HERO, the undeformed versions of these quantities are input by the user and 
for optimization problems represent the design variables. The blade layout routines 
require the deformed quantities which are obtained by superimposing the deformation 
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vector, s, upon the undeformed shape. Table 4-1 summarizes the information contained 
in the vector, fi. Denoting the variables at the left hand end of a segment by (*)l, those at 
the right hand end by (»)r, and the undeformed geometry by (»)o, then the deformed 
geometry, 


SL 

= SLo + U L - Ur 

SWEEP 

= SWEEPo + A V SO*cos(COLL) + AVGO*sin(COLL) 

TWG 

= TWG 0 + [(Rx)l-(Rx)r]/sl 

TWR 

= TWRo + (Rx)l 

ANH 

= ANHo + AVGO*cos(COLL) - AVSO*sin(COLL) 

CHORD 

- unchanged by s 

COLL 

- unchanged by s 


(4-40) 

where, AVSO = 

\ [ (Rz)l + (Rz)r] 

AVGO = - 

\ [ (Ry)l + (Ry)r] 


This takes into account the global rotation about the X-axis of the entire blade when 
applying the collective after the blade planform has been defined. 


4.6 A Note on the Wake/Structure Coupling Matrix 

The structural equations, Eq. (4-37) can be expressed as. 


Af={ f-K S ) = KAs (4-41) 


where the deformation is updated s+A$. 
equations so that. 
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Qqx Qqy ® 
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Qwx Qw y 0 
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Ay • 

Af 


o 
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As 


This can be combined with the wake 

(4-42) 


The left hand side residuals formed from nonlinear equations are driven to zero by an 
iterative process that amounts to a Newton method for updating the Xo y, and £. The right 
hand matrix of first order linear derivative terms does not reflect the coupling between the 
aerodynamic and structural equations and so the question arises as to whether the 
coupling terms should also be derived. Note that the accuracy of the analysis is entirely 
dependent upon the computation of the left hand vector of Eq. (4-42). The implicit 
coupling of the two sets of equations is also embodied in the left hand terms. The role of 
the right hand side matrix is to enhance the convergence of the iterative scheme providing 
superior stability and rate of convergence to the overall scheme. Thus, provided that this 
goal is achieved, one is at liberty to make approximations to its elements. 

Our experience indicates that neglecting the block off-diagonal coupling matrices is 
both an adequate and desirable approximation. Earlier version of the code included the 
off-diagonal coupling terms which were derived at considerable computational cost. It 
was found that the overall rate of convergence was minimal when compared against the 
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present case with matrix coupling entries neglected. Furthermore, the robustness of the 
overall solution scheme was comparable in both cases and was found to be best 
controlled by specifying the parameters DLD and DLR in sensitive cases. The 
approximation is very desirable for several reasons: 1) Evaluation of the off-diagonal 
terms is very expensive in computer time and requires substantial storage space; 2) The 
advantages of the efficient banded LU decomposition scheme described in Section 4.4 are 
forfeited. Instead, the Gaussian elimination scheme employed in EHPIC/HERO for 
solving the N wake equations is now applied to an (N+NF)*(N+NF) array where NF is 
the no. of degrees of freedom in the structural model. Hence the number of operations 
involved in Gaussian elimination increases from 0(N 3 ) to 0([N+NF] 3 ). 3) There is a 
substantial reduction in coding complexity. 4) For optimization studies, the effects of 
structural deformation are more efficiently incorporated by recomputing the structural 
stiffness properties and deflections for perturbed designs. Therefore, the coupling entries 
in the right hand matrix term in Eq. (4-42) may be safely neglected resulting in 
substantial reduction in computational resources with no noticeable loss in robustness. 
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TABLE 4-1 


Specification of the Global Degrees of Freedom For Element, ie. 


Global Deformation Degree of Freedom 


Left-hand node of element, ie: 

V 

s 8ie-7 


Rz 

s 8ie-6 


w 

s 8ie-5 


Ry 

s 8ie-4 


Rx 

s 8ie-3 


U 

s 8ie-2 

Mid-node of element, ie: 

qio 

s 8ie-l 


qi3 

s 8ie 

Right-hand node of element, ie: 

V 

s 8ie+l 


Rz 

s 8ie+2 


w 

s 8ie+3 


Ry 

s 8ie+4 


Rx 

s 8ie+5 


u 

s 8ie+6 
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5. IMPLEMENTATION OF DESIGN OPTIMIZATION 
5.1 Outline of the Optimization Solution Method 

The fundamental procedural elements involved in optimizing the rotor design in 
EHPIC/HERO are outlined in Fig. 5-1. The initialization routine specifies the objective 
function, the design variables and the optimization method to be employed. The 
subsequent optimization process is iterative and involves sequential evaluation of local 
objective function and constraint derivatives, pre-processing this information into suitable 
input to the optimization algorithm, implementing the improved design, and re- 
converging the wake analysis in EHPIC/HERO. The basic routines that carry out these 
procedures are also indicated in Fig. 5-1. This section describes each of these operations 
and their computational implementations. 


5.1.1 Optimization algorithms 

It is appropriate to begin with a description of the optimization algorithm since the 
input required to specify the minimization problem dictates the information that must be 
supplied by other routines. Thus, the set of computer code modifications and additions 
that have been effected in EHPIC/HERO to carry out the optimization process depend 
fundamentally on the core optimization algorithms employed. Conversely however, for 
problems of this scale, the computational cost involved in providing this information will 
limit our choice of practical optimization methods. Therefore, the selection of a suitable 
optimization scheme is itself an attempt to optimize the fundamental trade-off between the 
rate of convergence to an optimal design over a series of iteration steps and the amount of 
computation required per iteration. 

In general, the constrained optimization problem may be posed: 

Minimize J(X) subject to the constraints 

gj(X) £ 0 , j=l,...,mj 

gk(X) = 0 , k=mj+l,. . .,m^ 

lj — Xj ^ Uj , i=l,...,n 

(5-1) 


where X is the state vector of order n, J is the objective function to be minimized by 
appropriate choice of X. the gj and g^ are the mj inequality and mg equality constraints, 
and lj and Uj define lower and upper bounds for the corresponding design variable, Xj. In 
the context of rotor design optimization in hover, the complete state variable vector, 

X = { 2Sc» d } T where Xc and y are the wake collocation point locations and blade 
quadrilateral bound circulations used in the EHPIC/HERO analysis, £, are the structural 
deflections, and the vector, d, constitutes the set of design variables such as blade twist, 
sweep, and anhedral. 

In practice, for the performance objective functions considered here, structural 
deformations, & are eliminated prior to posing the optimization task. Instead, effects due 
to structural deformation are implicitly accounted for when evaluating the derivatives for 
the design vector, d- Candidate objective functions, JQQ, include the power arising from 
induced drag or profile drag, thrust, the Figure of Merit, or a combination of these. 
Examples of imposed constraints include such requirements as maintaining constant 
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thrust or remaining below a specified total power ceiling. Note that the equations used in 
the wake calculation are used to form equality constraints, gk, relating the 2c> X and d. 
[see Eq. (5-26)]. As described below, these equality constraints can be either used 

directly as input into the optimization routine or else employed to eliminate the 2c, X The 
latter is more efficient, but one forfeits the ability to impose constraints on the 

perturbations in 2c and y. 


In general, the functions, J and gi are nonlinear. Furthermore, analytical expressions 
for these functions in terms of X are usually unavailable and their evaluation is carried 
out in a post-processing phase after the hover solution has converged. Thus, although 
J(X) and g(X) can be obtained given X the derivatives of these functions are typically 
available only by numerical differencing. We observe that the influence coefficients 
employed in the wake analysis are basically finite difference approximations to the first 
derivatives of collocation cross-flow velocities, q, and downwash velocities, s', wj.t. 2c, 


and y. The approximation of the first derivatives of these quantities wj.t. to the design 
parameters, & exacts a computational cost comparable to that of obtaining the other 
influence coefficients and is judged acceptable. Thus, an optimization technique that uses 
first order information is highly desirable since much of the necessary data is already 
available from the wake analysis, and the remaining information can be determined with 
the same order of computational cost as incurred in one wake analysis iteration. Zeroth 
order techniques using only the evaluated J(X) and gi(X) are wasteful of the first order 
available information and require excessive iterations to converge. Second order methods 
are judged prohibitive as a result of excessive computational cost involved in evaluating 


the second order derivatives, 


d 2 J 



and 



IA=A 


5.1.2 Review of the Phase I optimization scheme 

The first version of the wake analysis with optimization capability employed a 
sequential linear programming (SLP) technique as the core optimization technique 
(References 51, 52). As described in Reference 4, the nonlinear, J(X) and g(X) arc 
linearized about the current design and wake solution, 

J(X) » J(X* ) + VJ(X* ) AX (5-2) 

where the constraint 

0 = g k (X)“gk(X*) + Vg k (X*)AX (5-3) 

is imposed. In Reference 4, JfX) was equated with the total power required by the blade, 
and the set of equality constraints, gi(X)» is composed of a thrust constraint to maintain 
constant thrust The set of equations used in the wake solution and extended to account 
for design variable perturbations. The twist gradient over each blade segment constituted 
the design vector. The state vector, X*, denotes the current design and converged wake 

description, and AX perturbations about that solution. Specifically, the posed problem 
was to minimize the power expressed as. 
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(5-4) 


P-P(X’) + {Px.Prjjly} + { p dK A d} 

subject to constant thrust, 

0 = T-T 0 - T(X*)-T 0 + {T I ,T r }|^j + {T d }{Ad} (5-5) 

the wake constraints. 



and move limits, lj <, Axj ^ uj , where, 


all 

II 

p? 

; P- 

_3P . p,. 
r '9* ' Pd ’ 

ap 

(5-7 a) 

X 37 

t * = 3£ 

; T 

dT „ 

; Td 

a t 

(5-7b) 

Qqd = 

1 1 

* Qwd = 

aw; 

adj 

(5-7c) 


and the entries of submatrices Qqx, Qqy, Q wx , and, Qwy are identically the influence 
coefficients employed in the EHPIC/HERO analysis. Thus the additional computational 
load involved computation of derivatives w.r.t to the design and the thrust and power 
derivatives. The wake constraints ensure that, to a linear approximation, the cross-flow 
and downwash velocities remain zero under the design perturbations. Eqs. (5-1) were 
solved using a simplex linear programming routine (Reference 52). The state is then 

updated using the computed AX and the EHPIC/HERO analysis reconverged. The 
process continues sequentially until no substantial improvement in JOO is obtained. 

The original SLP algorithm made efficient use of the first order information at hand 
and proved sufficient in improving the blade twist distribution for power reduction at 
constant thrust. Nonetheless, certain improvements to the algorithm were desirable. 
Chief among these were: a more efficient implementation of the upper and lower limit 
constraints which originally were entered as 2n inequality constraints (n being the order 
of the state vector, X); a means of converging to minima that lie off the imposed 
constraints; improved numerical conditioning of the simplex tableau; and ways of 
reducing the tableau array size. Hence, an integral task of this effort has been evaluating 
alternate optimization algorithms and addressing ways to improve the efficiency of the 
algorithm. 


5.1.3 Selection of optimization algorithms for Phase II 

As explained above, our examination of optimization algorithms has been limited to 
those requiring at most first order information. The available techniques may be broadly 
categorized into constrained and unconstrained minimization methods. The former 
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approach deals with the imposed constraints directly, whereas the latter augments the 
objective function with a penalty function which basically penalizes violation of the 
imposed constraints. More stringent constraints can be emphasized by varying the 
weighting parameters in the penalty function. The transformed problem is then solved 
using an unconstrained optimization method. Unfortunately, the unconstrained 
techniques tend to suffer from numerical ill-conditioning related to the penalty function 
weightings. Furthermore, it is difficult to judge a priori which of the constraints are of 
greater importance in the problem and how to quantify the weighting terms to 
appropriately enforce a given constraint Finally, computational experience in structural 
optimization problems indicates that such schemes frequently converge much more 
slowly when compared to the direct methods (References 51, 53). 

Therefore, we focused on the constrained optimization methods. The candidates 
considered were an extension of SLP with move limit reduction near the optimum, the 
method of feasible directions, and sequential quadratic programming (SQP) (Reference 
51). It was decided to implement SQP for the following reasons: 


(1) It tends to exhibit superior convergence properties. This is due to the iterative 
generation of a Hessian matrix, which is essentially a positive definite approximation to 

^2j 

. As detailed below, the matrix can be generated using only first order 

x=x* 


dX 1 


derivatives. Qualitatively, it embodies information from previous steps to estimate the 
curvature, or second order properties, in the vicinity of the optimum. It is this 'memory' 
of previous iterations which gives SQP methods an advantage over other techniques 
(SLP, method of feasible directions) using only the local first order data and which 
originally engendered a closer look at SQP methods. 


(2) The presence of a positive definite quadratic term in the objective function enhances 
convergence near minima that do not lie on the constraint boundaries. It is well known 
that the solution to a linear programming problem lies at the intersection of exactly n non- 
degenerate constraints. Thus, if the true nonlinear minimum does not lie at such an 
intersection the linear programming solutions will dither between the imposed move 
limits, li and ui. In practice, the move limits are reduced when proximity to the minimum 
is detected so that one effectively 'shrinks' the feasible domain down to the optimum 
design. This procedure is time consuming and is circumvented in SQP methods by virtue 
of using a quadratic objective function. 


(3) The SQP problem is solved by using a modified version of the simplex algorithm. As 
described below, solving a quadratic programming problem amounts to finding a feasible 
solution to the Kuhn-Tucker conditions which can be cast as a linear programming 
problem with additional logical constraints in the selection of variables leaving and 
entering the active set (i.e. the set of variables that away from its bounds, lj and ui.) Thus, 
the core simplex routines employed by SLP and SQP methods share a strong commonalty 
and in fact, the same simplex algorithm is used by both methods in the optimization code 
developed here. 


One drawback of the SQP approach is increased array dimension of the simplex 
tableau. Whereas the tableau for the modified SLP scheme developed here is of order (m 
by n) where m is the number of imposed constraints, the SQP solution process requires an 
array of order (m+n by m+2n). Since the number of wake related equality constraints is 
n, then m=n and the SQP method has six times the memory requirements of the SLP 
method used to solve the same problem. The number of floating point operations 
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increases accordingly. Furthermore, away from the optimum, the SLP and SQP 
algorithms tend to advance at the same rate. This is due to the fact that away from the 
optimum the design change per optimization step is limited primarily by the imposed 
move limits. The linear approximation to the objective function dominates the quadratic 
terms and the design changes obtained per iteration tend to be identical using either SLP 
or SQP. The picture is different in the neighborhood of the minimum where many of the 
first order derivatives tend to zero, so that the quadratic term in the SQP formulation 
becomes important and accelerates convergence to the final design. Accordingly, it has 
been found to be desirable to retain both the SLP and SQP techniques, using SLP initially 
and then continuing with SQP when proximity to the optimum is ascertained. The 
similarity of the underlying simplex procedure used in both SLP and SQP solution 
procedures has thus proven to be a distinct advantage. 

In the following sections, the linear programming (LP) and quadratic programming 
(QP) problems are posed, and the solution to the QP problem stated and its relation to the 
LP problem established. The construction of quadratic term in the SQP approach is then 
described and finally some aspects of the common simplex routine detailed. 


5.2 Statement of the SLP and SQP Problems 

The nonlinear optimization problem posed in Eqs. (5-1) is repeated: 

Minimize J(X) subject to the constraints 

gj(X) ^ 0 , j=l,...,mj 

gk(X) = 0 , k=mj+l,...,m e 

lj^xj^Uj , i=l,...,n 

(5-8) 

The solution at iteration step, q, is known, and the problem is to determine a change 
in this solution which results in a reduction in J(X). In the iteration limit, convergence to 
a local minimum is expected. The LP problem to be solved each iteration step is: 

Minimize, J(X q ) + VJ(X q ) AX subject to 


gj(X q ) + V gj (X<l)AX<0 



g k (X q ) + Vg k (X q ) AX = 0 

, k = iD|‘^ 1 , . . • 


Alj < Axj < Auj 

» 1— 

(5-9) 

The upper and lower bounds incorporate move limit requirements: 


Ali = Max. { - Ax™“ , lj - x? } 


(5- 10a) 

Aui = Min. { Ax™** , Ui - xf } 


(5- 10b) 


where Axj”** is the maximum change in xi allowed per iteration step. The QP problem is 

identical to the LP problem except that the objective function is augmented by a quadratic 
term to: 
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(5-11) 


J(X q ) + VJ(X q )AX + |ax t [b] AX 


where matrix, [B], is a positive definite approximation to the Hessian of the objective 
function and local constraints and is generated iteratively from successive evaluations of 
the first order derivatives as described in Section 5.2.2. 

The LP problem at each iteration of the SLP method is solved using the simplex 
technique described in Appendix A. In the SQP case, the sequence of QP problems is 
solved finding feasible solutions to the Kuhn-Tucker equations as described below. 


5.2.1 Kuhn-Tucker conditions 

The solution to the QP problem is stated in terms of the well-known Kuhn-Tucker 
conditions (Reference 5 1) as follows. The full nonlinear conditions for an extremum are 


thereby stated by first defining a Lagrangian, 



n n 

mi+me 


LQD = J(X) + 5>,j(lj - Xi) + 5>2i(Xi - Ui) + 

2 

Aj &(X) 

i i 

i=l 

(5-12) 

then in addition the imposed constraints. 



VL(X*)=^r = 0 

dx " 


(5- 13a) 

V\\ (if ~li) = 0 ; Hx (ui -X-) = 0 ; jiJ , /4 

£ 0 

(5- 13b) 

Aj gj(x*) = 0 ; X* £ 0 (j=l,...,mj) 


(5- 13c) 

A k unbounded (k=mj+l,...,mi+mc) 


(5- 13d) 


where the superscript, (•)*, denotes evaluation at the optimum point. The new 

parameters, jii, Ji2 and, are Lagrange multipliers that correspond to the lower and upper 
bound constraints and the imposed constraint equations respectively. The first of these 
equations define an extremum while the remaining conditions restrict the allowable 
values for the Lagrange multipliers. Specifically, for inequality constraints, the 
associated Lagrange multiplier can only be non-zero if the constraint is active. This 
implicitly states that when a given inequality constraint is active then the gradients of the 
cost function and constraint point in opposite directions so that the design cannot be 
further improved without violating the local constraints. 

For the SQP problem posed in terms of design perturbations, in addition to satisfying 
the imposed constraints, the optimum AX satisfies: 
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with 


mi+me 

VJ(X q ) + [B] AX - M t + M 2 + X ^i v «i(x q ) = Q (5- 14a) 

i=l 


A*li (Axi — Alj) = 0 ; /i 2 i (Auj — Axi) = 0 ; 

^li • /*2i ^ 0 

(i=l,...,i 
(5- 14b) 

*j[gj(x q ) + V gj (x q ) Axj = 0 ; Aj 2 ! 0 

f 

* 

II 

(5- 14c) 

Xk unbounded (k=mj+ 1 , . . . 411 c) 


(5-14d) 


In the above, the first n equations are merely linear equality constraints which are directly 
accommodated in the simplex formulation. Likewise, the positivity requirements upon 

the components of iii , and. Xj are naturally dealt with in the LP framework. The 
imposed linearized constraints, gj, are of course already handled by the simplex 
algorithm. Therefore, the only modifications required to solve the QP problem utilizing 
the simplex algorithm are those related to the nonlinear constraints in Eqs. (5-14b,c). 
Fortunately, these conditions can be implemented entirely by logical restrictions upon the 
variables allowed to enter and leave the basic set during a pivot operation. Specifically, 
in the simplex formulation the inequality constraints, are converted to equality constraints 
by the use of slack variables, Sj: 

gj(x q ) + Vg j (x‘l)AX £ 0 (5- 15a) 

gj(x q ) + Vg j (x‘l)AX - sj = 0 , sj £ 0 (5- 15b) 

Thus the related Lagrange variable, Xj, can only enter the basic set if Sj is at its bound, 
Sj=0, implying that the corresponding constraint is active. Equivalently, only one of sj 

and Xj can be in the basic set at a given time. Since membership of a given variable in the 
basic set is a logical attribute, the nonlinear constraints embodied in the Kuhn-Tucker 
conditions are implemented simply by conducting a series of membership tests upon the 
candidate pivot elements. 

The problem solved by the QP algorithm is an approximation to the original nonlinear 
optimization task posed in Eqs. (5-8). Thus in the spirit of iterative methods for solving 

nonlinear problems by local linearization, the solution, AX. obtained is used to update the 
current design, the functional and constraints and their first order derivatives re-evaluated 
about the new state, and a new QP problem formulated about that point. Although it is 
difficult to derive conditions for convergence in the general case, the sequential process is 
expected to converge to the nonlinear optimum, especially if the functional is convex in 
the vicinity of the optimum. 

An extension to the SQP procedure described is to employ the state vector update, 
AX, as a search direction and perform a 1-D minimization or equivalently determining an 
optimal step length along this direction. The theoretical advantages of conducting a 1-D 
search along the optimal direction include more accurate update of the Hessian matrix 
approximation. Also, theoretical proofs of superlinear convergence of the SQP scheme 
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assume that 1-D optimization is executed. Unfortunately, at least two additional 
functional evaluations (calls to EHPIC/HERO), are necessary to determine an 
approximate 1-D minimum using a quadratic polynomial fit. This cost was found to be 
unnecessarily expensive since the rate of convergence to an optimum was not 
significantly improved. We believe this to be due to violation of constraints that can 
occur when the optimal step length differs from unity. It was also observed that the 1-D 
functional which is constructed from the cost function and active constraints is sensitive 
to round-off error and constraint violation and occasionally leads to unreasonable step 
lengths which must be limited. Furthermore, as noted earlier, when the design change is 
limited primarily by the imposed move-limits, the progression of the solution toward the 
optimum is fairly constant whether employing SLP, SQP or SQP with 1-D minimization. 
Hence, the SQP algorithm employed here is executed using unity step length. 


5.2.2 Hessian matrix update formula 

The Hessian matrix approximation, B, is updated at each optimization step, q, 
according to the Broydon-Fletcher-Shanno-Goldfarb update formula (References 51, 54, 
55). Given the current approximation, BQ, the updating proceeds as follows: First rescale 
BQ as. 


B* = 


k = 


p T w q 

P T B q P 


then update, 


where. 


B „l = B - _ bWb' + 

PB p 

P s X q -X^ 1 = AX q_1 



r = 0w q + (1 - 0)B*p 


w q = V x L(x q ,i q -') - 

L = J(X) + £ Aj g|(X) 
iel h 


and 


1 


e = 


0.8 p T B*p 

Ip b p - p 1 t 


if p T r ^ 0.2p T B*p 
if p T r < 0.2p T B*p 


(5-16) 


(5- 17a) 

(5- 17b) 
(5-17c) 
(5-17d) 
(5-17e) 


(5-17f) 


Here, Ih is the set of holonomic constraints as explained in the next sub-section. The 

matrix, B q , is first scaled by the factor k as recommended by Luenberger (Reference 56). 
For a quadratic prog rammin g problem where the true Hessian of the Lagrangian function 

at the optimum is H, then defining k in this manner guarantees that the range of 
eigenvalues of BH' 1 spans unity at each iteration. Furthermore, the condition number of 
(BH" 1 ) will be non-decreasing which is desirable from a numerical standpoint. In theory, 
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the specification of q above ensures that if B is initially positive definite then it remains 
so after the updating procedure. However, this is not necessarily true in the presence of 
round-off error. Specifically, the last term in Eq. (5- 17a) is always positive semi- 
definite. However, numerical differencing between the first two terms can introduce 
sufficient error to render B*!' 4 ' 1 non-positive definite. Hence, the expression is modified 
slightly by adding small positive terms to the diagonal elements of B: 

B q+1 <- B q+1 + Tg*p [diag {y?}] . y = B *P (5-18) 

The entries of this additive diagonal matrix are identical to the diagonal components of 
the second term on the R.H.S. of Eq. (5- 17a) scaled by a factor of 0.01. The manner in 
which information of the cost function and constraint derivatives from preceding 
optimization steps is implicitly incorporated into BQ is apparent from the definition of w4. 
During the course of the optimization, round-off error and higher order contributions 
from die constraints and cost function derivatives are also accumulated in B. Thus, 
periodic resetting of the B to a scaled identity matrix is recommended and is done in the 
code every 2n-2 optimization steps (resetting B every n steps can occasionally slow 
convergence) and also when convergence to the optimum is detected. 

The preceding construct for the matrix B guarantees that the generated matrices B9 
remain positive definite throughout the optimization provided that the initial matrix, 
B°>0. This is true regardless of whether the true nonlinear cost function, J, contains 
maxima or saddle points. Under certain assumptions on the convexity of the objective 
function and constraints (References 54, 55) it is proven that supcrlinear convergence 
results when using this update formula. The sequence, BQ, need not converge to the 
actual Hessian of the Lagrangian at the solution. Instead, the projection of B<? unto the 
tangent vector space generated by the linearized acting constraints converges to the 
corresponding projection of the true Hessian. Or symbolically, if the optimum solution is 

denoted by ( X*J1*A* ) then 

lim Z t [v xx l(x*,/£*, A* j - B q ]z = 0 (5-19) 

where, the Lagrangian is defined in Eq. (5-12), and Z is any vector lying in the space 
generated by the active constraints: 

VgjQD , jejA and Vg^QQ , k=l,...,me (5-20) 

where the set of active inequality constraints, jA={j: gj=0 , j=l,...,mi). 


5.2.3 Distinction between holonomic and non-holonomic constraints 

The set of constraints, Ih, appearing in the summation of the constraints in Eq. (5- 
17e), is the subset of the complete set of constraints that are holonomic. The update 
formula for B9 assumes constraints of the form Eq. (5-8) which are functions of state, X, 
and possibly other parameters which are fixed during the optimization process. These 
constraints may be thought of as global constraints in that they form well defined and 
fixed hypersurfaces in the state space. Constraints may also be stated in differential form 


66 



which in the present context implies that the constraints arc expressed in the form, dgj(X, 
A£)<Q or dgkQL AX)=0 which may be regarded as local constraints imposed at a given 
optimization step but which may change between steps. When such constraints arc 
integrable so that they can be expressed in the form gjGD^O and gk(2D=0 respectively 
then they are termed holonomic. When no such integrated form exists, they are classified 
as non-holonomic. 

Examples of holonomic constraints applied to the hover optimization task include the 
specified upper and lower bounds on 2L limits on the nodal twist, or equality and 
inequality constraints on power and/or thrust Examples of non-holonomic constraints 
include move-limit constraints, changes in nodal twist and, most significantly here, the 
equality constraints derived from the wake analysis. The wake equations are expressed in 

terms of the changes in collocation point positions, A&c and bound circulations, Ay. 
These perturbations are determined in the optimization analysis and used to update 2c and 

y. However, during the subsequent EHPIC/HERO analysis where the wake is re- 

converged, 2 c and y will generally change, i.e., the values of 2c and y differ between calls 
to the optimization process with the discrepancy being due to nonlinearities. Therefore, 
these constraints are not included in the update formulas. Similarly, the move limit 
related inequality constraints are not used for the Hessian matrix update. However, when 
the solution is sufficiently close to a 'hard' bound (upper or lower limit imposed on the 
variable) so that it forms one of the boundaries of the feasible region then the associated 
inequality constraint now is used in the update equation for B4. Naturally, all of the 
constraints must be retained in the SQP optimization in order to define the feasible 
region. 


5.2.4 Extensions to SLP 

The SLP algorithm utilizes the same simplex subroutines as the SQP scheme. 
Beyond the modifications made to the simplex algorithm itself, the chief modification 
that has been made from the Phase I algorithm is additional logic for the sizing of the 

move limit, Axj 112 *, associated with each xj. The rationale is to retain constant size move 
limits as long as the design point is heading towards the optimum. If the optimum 
solution is fully constrained (i.e. lies at the intersection of a total of n linearly 
independent constraints) all is well since the LP solution will then coincide with the 
actual nonlinear solution. In general however, the optimum solution may not be fully 
constrained and when one is sufficiently near the optimum such that the optimal design 
point lies within the feasible region, the sequence of points obtained via SLP from there 
on oscillates about the optimum essentially bouncing between the move limits defined for 
the design variables. Thus one has obtained the optimum solution only to within the 

resolution defined by the Ax f 13 * and to increase the accuracy of this solution the 

Ax I" 8 * must then be reduced. This is easily done if one is able to ascertain when the 
optimum lies within the feasible region surrounding the current design point. This is 
accomplished by keeping track of the design changes. Specifically, the two most recent 
design points are stored and if all of the design variables are observed to remain within 

the region defined by X ± A2L max then the move limits are reduced, i.e., if 
IxJ-x^Axf 18 * , Vi, then set Ax}" 8 * <- 0.55 Axj 11 **. 
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5.3 Application to Rotor Design Optimization 

The implementation of the sequential optimization techniques outlined above, 
amounts to specifying a design vector, d, together with the desired cost function, J, 
constraints, gi, and their derivatives as input to the optimization algorithm. The design 
vector, d, is determined by the user in the optimization input file and consist of the 
parameters that specify the undeformed blade geometry. These presently include sweep 
distribution, twist distribution, anhedral distribution, chord distribution, and collective. 
The coding for blade radius and cutout optimization is also in place but their use as 
design parameters is not recommended at present. The equations used in the 
EHPIC/HERO wake analysis form a set of equality constraints in the optimization 
analysis. In many cases, the equality constraints associated with the wake equations can 
be solved directly to eliminate some of the variables thereby greatly reducing the 
dimensions of the optimization tableau. These topics together with the objective 
functions and constraints that may be applied are examined in the following sections. 


5.3.1 Description of available objective functions 

The user is offered a selection of objective functions to be minimized. In the present 
version these are restricted to thrust, induced power, profile power, or combinations of 
these. Objective functions associated with structural properties can also be formulated to 
mi nimi ze ma ximum stresses or deflections. However, structural optimization has not 
been included as an option for several reasons. First, the focus of this effort has been 
upon the aerodynamic aspects of the hover problem and the purpose of including a 
structural modeling capability is to improve hover prediction for real rotors. To that end, 
a structural model that predicts the blade deflections with an accuracy that is consistent 
with the amount of detail contained in the geometry specifications, is adequate. Hence, 
for example, the number of finite elements is commensurate with the number of blade 
segments. This also allows a concise structural specification file where the major blade 
stiffness and mass properties are characterized in terms of area integrals. Finally, the 
model can be constructed very efficiently which is advantageous when determining the 
design related derivatives in the optimization analysis since the stiffness matrix and nodal 
forces must be recomputed for each design perturbation in order to account for 
deformation effects (see Chapter 4). Second, it is judged that a serious approach to the 
task of structural optimization requires a more detailed structural model in order to obtain 
accurate estimates of the blade stresses. This requires specification of skin thickness, and 
other localized geometry information rather than the integrated cross-section parameters 
used in EHPIC/HERO. Finally, the structural optimization task is constrained by other 
considerations not directly relevant to the hover problem. Notably, the modal 
frequencies and blade dynamic response in forward flight are expected to be more 
significant in guiding structural optimization. 

Therefore, we have sought to furnish a series of objective functions which represent 
integrated aerodynamic performance parameters deemed to be of most interest to 
potential users. At present the fundamental parameters are the thrust, T, induced power, 
Pi, and profile power, P p . The user selects the objective function by specifying the input 
parameter KOBJ(l) as follows (see Reference 12): 

KOB J(l) Operation Performed 

1.11 Minimize Total Power, Pt = Pp+ Pi 

2. 12 Maximize Thrust, T 
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3 Maximize Figure of Merit, V2~ Cy 3 / 2 / 2 Cq. 

4 Maximize Propulsive Efficiency, Ct/Cq. 

20 Multi-objective minimization of both thrust and total 

power 

30 Multi-objective optimization of combination of thrust 

and induced and profile power contributions. 


where Cr=T/(p7i£2 2 R 4 ) and CQ=(PffPp)/(pfl£2 2 R 5 ). The cost function combinations for 
KOBJ(1)>20 are specified below. It suffices here to note that these involve combinations 
of T, Pi, and Pn. Then it follows that all options for KOBJ(l) are constructed from one or 
more of the functions, T, Pi, and Pj>, which must be supplied in the EHPIC/HERO 
analysis along with the first order derivatives which are calculated by finite differencing. 

Except for options KOBJ(l)=3,4, the objective functions appear to exhibit weak 
curvature so that the design tends to progress indefinitely along a path of steepest descent. 
This is similar to the problem encountered in linear programming procedures where the 
objective function is linear and unless suitably constrained, can return optimal designs 
with infinite values. It is clear, for example in a power optimization calculation, that 
unless a lower bound upon thrust is imposed, the design will move in a direction that 
reduced thrust along with the power leading to zero or negative thrust levels. To properly 
pose the optimizations, each option of KOBJ(l) also implicitly adds constraints: 


KOBim 


Applied Constraint(s) 


1 

2 

3,4 

11 

12 

20 

30 


T - Tgnec 
Py = (Py)spec 
None 
T^T 


spec 


Pt ^ (PlOspec 
T ^ Tgpec and Py ^ (Pr)spec 


T ^ Tjpec » Pi ^ (Pi) spec » and Pp ^ (Pp)spec 
where (*)spec is a specified level (default levels correspond to the initial design). 


5.3.2 Choice of multi-objective functions 

A multi -objective design task is generally approached by constructing a global 
objective function from the single objective criteria. The simplest form is simply a linear 
combination of the constituent cost functions and utilize this as the cost function, J, in the 
same optimization procedure employed in the single objective case. A preferable 
approach is to combine the single objective functions in a possibly nonlinear fashion so 
that the resulting global cost function represents a physically meaningful entity. 
Examples of this are the options KOBJ(l)=3, 4 where the thrust and power are combined 
nonlinearly to form the figure of merit and propulsive efficiency. Often, however, such 
convenient combinations are not obvious, particularly in multi-disciplinary optimization. 
Therefore, unless or until global optimization cost parameters are available, the linear 
combination of single objective functions is arguably most useful due to ease of 
implementation and tractability, and since sensitivity analyses which help gauge relative 
trade-offs between the various constituent cost can be straightforwardly conducted for 
this combination. 
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If the single objective cost functions are denoted by Ji then the composite cost 
function. 


J G 


X w i J i 


(5-21) 


where the wj are weighting parameters that are intended to reflect the degree of 
importance attached to minimizing the constituent, Jj. The specific forms for the multi- 
objective options of KOBJ(l) are: 


KOBJ(1)=20: 
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where T= thrust, Px=total power, Pj=induced power and P p =profile power. The (*)mi n 
and (*)max values are either specified by the user or set to default values in 
EHPIC/HERO. Their purpose is to scale the dimensional quantities so that the Jj have 
roughly equal influence upon the overall Jq. Ideally, if the wi were all equal, then the Ji 
would also be equal in the globally optimal design. In practice best guesses to these 
ranges are usually made. One alternative (Reference 57) is to perform the single 
objective function minimizations separately. For each design thus obtained, dk say, 
EHPIC/HERO may be used to evaluate T(dk), Pj(dk), and P p (&k). Then 

T^ = max{T(d k )}, and similarly for T m in, (Pj)max etc. 

An important consideration is the selection of the wj. This is a non-trivial task with 
the goal of choosing values for wj that leads to an overall qualitatively 'best' design. It 
often happens that one of the Ji dominates and the optimal design is similar to the 
corresponding single objective case. A more serious effect is that usually some of the Jj 
are 'sacrificed* in order to allow decrease in some other Jj so that a net decrease in Jg is 
attained. Problems arise if the change in the numerical magnitude of Jj is small in 
comparison to the other Jj, but that this change represents an unacceptable increase of that 
cost. Essentially, this reiterates the desirability of a global cost function which is derived 
from physical considerations (e.g. figure of merit; propulsive efficiency) since it is 
otherwise not clear what Jg is intended to represent Hence it is assumed that this issue 
has been carefully considered by the user who is required to specify the Wj in the input 
file. 

One additional safeguard is supplied however which prevents the individual costs 
from exceeding specified upper bounds. Thus the following inequality constraints are 
enforced for, 

KOBJ(1)=20: 
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In fuzzy optimization terminology, the Ji may now be regarded as membership functions. 
The combination of the normalization terms and the limiting effected in the preceding 

equations implies that the Jj € [0,1]. It is assumed here that by suitable choice of (*)max 
and (*) m jn the individual Jj^O. Thus the individual cost functions are placed on an equal 
footing with the other imposed constraints. This guards against the 'sacrifice' behavior 
mentioned above and tends to produce qualitatively better designs that give overall 
improvement in at least some of the Jj while ensuring that the others do not become 
unacceptably large. 


5.3.3 Incorporation of the wake constraints 

The equality constraints upon linearized normal/binormal velocities and the blade 
quadrilateral downwash velocities are directly available from the EHPIC/HERO analysis 
and are expressed as, 

fAqjfOl Qqx Qqy Qqd 

(AwJ [Oj Qwx Qwy Qwd 

These constraints play two roles in the optimization analysis. First, when enforced, they 
guarantee that to within a linear approximation, the post-optimization design maintains 
zero downwash and collocation cross-velocities as required for the converged wake 
solution. In practice, the inherent nonlinearities require that the wake solution be 
reconverged and the post optimization design solution provides a very good first guess to 
the converged state. Second, these constraints serve to linearly couple the design 
perturbation with the wake changes so that the total influence of a given design variable 
upon the constraints and cost is accounted for. In EHPIC/HERO, the design related 
influence coefficients are effectively partial derivatives with the collocation points and 
blade circulations held fixed. Thus, for example, the thrust coefficient, T<j, for design 
variable, dj, is computed by perturbing the blade geometry corresponding to dj and then 
forming the finite difference of the thrust Thus, 

Td = gjl Ok, Affixed) (5-27) 



(5-26) 
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To determine the physical change in thrust that occurs when d is perturbed, the total 
derivative must be evaluated implying that the changes in the wake, Axe, and bound 
circulations, Ax must be taken into account. Thus, 

dT dT . ^dT . dT 

dd = 3*** 3 ^ 5T ^ (5 ’ 28a) 

or, 

g = T, tec + TyA* + T d Ad (5-28b) 

The wake constraints can now be used to solve for the wake position and bound 
circulation perturbations that take place for perturbation Ad* 


to l-J 
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Qwx 
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_Qwd_ 


so that the physical thrust change due to Ad is, 


AT = 


Td-fr, T r } 


Qqx Qqy 

pi 


Qwx Qwy 


[Qwd j ^ 


Ad 


(5-29) 


(5-30) 


This is extremely important since in general the second term in (•) of the preceding 
expression dominates implying that the chief contributions are due to the related wake 
and bound circulation changes. In fact, it is the matrix partition, Q W g> which accounts for 

most of the contribution to AT. This has been observed to be generally true of 
aerodynamic parameters other than thrust (e.g., power). 

The set of linear equality constraints, Eq. (5-26), can either be entered into the 
optimization analysis directly in forming the constraint set, g, or else used to eliminate 

the variables As c and Ax via Eq. (5-29) in a pre-processing step prior to the optimization 
routine. The advantage of the latter approach is a drastic reduction in the dimension of 
the simplex tableau since when Eq. (5-26) is entered directly in the optimization routine, 
it consumes most of the tableau. The drawback is that move limits can no longer be 

imposed upon Axe and Ay and so the move limits imposed upon the design variables may 
need to be reduced to ensure that the subsequent wake analysis reconverges. Reduced 
move limits imply slower progress toward the final design. 

The EHPIC/HERO optimization parameter, KOPT, allows the use to select between 
these options: 

KOPT=l - Use the full wake equations to define equality constraints in the 
optimization routine; 

KOPT=2 - Use the full wake equations to eliminate Ajfc and Ax and update the 
design variable related influence coefficients (e.g., Td, etc.) 
KOPT=3 - Same as KOPT=l except that only the downwash equations are 
retained and Q wx is assumed zero. Thus the constraints 
corresponding to the wake analysis reduce to: 
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K0PT=4 - Use the reduced equations of KOPT=3 to solve for Ay. Then, 

Axe = Q 

and Ax = - [ Qwy ] t Qwd 1 Ad 

The particular option to use will depend on the specific blade configuration and wake 
analysis used. Increasing KOPT results in fewer computations but tends to be less robust, 
particularly for analysis involving a large number of wake points. 


5.3.4 Incorporation of structural deformations 

The constraints equations, Eqs. (5-26), and the objective function definitions assume 

that those influence coefficients identified with the design perturbations, Ad, take into 
account any structural deformation effects. When structural deformation is present, one 
distinguishes between the undeformed geometry which is defined by the design variables, 
and the deformed geometry used in the wake analysis. 11115 is most important when 
determining the various derivatives with respect to the design variables (e.g., Qwd. id. 
etc.) by the finite differencing technique. If we denote here by do any variable used to 
specify the undeformed geometry and by ds. its value in the deformed case, then is it not 
generally true that perturbing do leads to an identical perturbation in ds- i e -. m general, 

Ado*Ads due to change and reorientation of the deformation vector. Furthermore, a 
perturbation in a single component of do alters the deformed geometry of the complete 
blade. Take, for example, the undeformed sweep of a given segment lseg, 
UNDSWP(iseg). Perturbing UNDSWP will reorient the force due to blade rotation, alter 
the stiffness matrix which has major contributions from rotational stiffening, and 
consequently change the deformation state and the deformed geometry for the whole 
blade Note also, that even if the actual aerodynamic loading were to remain constant the 
reorientation of the blade results in a different nodal force vector and hence deformation 

state. The difference, Ado-Ads. naturally depends upon the material blade stiffness. 
However, this difference can lead to significant changes in some of the derivatives w.r.t. 
to design, (*)d. This is especially true of the twist gradient, where the tip twist change can 
be twice that due to the rigid body perturbation alone for typical blades, thus causing 
readjustments in these (*)d on the same order as ( # )d itself. This corresponds to the 
classical theory of torsional divergence in aeroelasticity. 

In principle, the state vector can be expanded to include changes in the deformation 
state, so that the equality constraints are augmented to: 


Aq 
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► = •< 
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Here, the residual Af = f - Ks, and the new influence matrix partitions, Q qs , Qws» Qfx> 
Qfy, and Qfd represent the change in the quantity corresponding to the first subscript due 
to perturbations in the state associated with the second. This approach was originally 
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implemented during this effort. However, as mentioned in Section 4.6, the calculation of 
the new influence matrix partitions is rather lengthy and cumbersome and was eventually 
dropped in favor of the approach below. 

Rather than conduct a linearized analysis involving the structural deformation 
changes, it proves far more efficient to perturb the undeformed design component, 
recompute the structural properties and loading, update the deformed geomeuy, and then 
calculate the various aerodynamic derivatives (Qqd, Td, etc.) based upon this geometry. 
The computation required to conduct the complete structural analysis and update the 
deformed geometry is a small fraction of that entailed in computing the wake related 
influence coefficients. Thus, in essence we define: 



(m) 

dd , 


(5-32) 


so that 2£c and % are held constant, but s is allowed to vary. Fortunately, this is a situation 
where the actual nonlinear analysis proves more efficient than the linear approximation 
without loss of robustness in the computation. 

In order to calculate the new deformation, the aerodynamic loading must be known, 
but this is not available until the deformed geometry is available. In theory, a time 
consuming series of sub-iterations is required where the deformation and loading are 
updated until convergence is attained. In the present implementation, the deformation is 
computed only once for a given design perturbation using the aerodynamic forces of the 
unperturbed design. This is justified in light of the fact that most of the change in 
structural deformation will be due to reorientation of the segment axes and corresponding 
changes in the local to global transformations matrices used in obtaining the stiffness 
matrix and nodal forces. Furthermore, this is consistent with the linear approximations 
used elsewhere in deriving the first order derivatives. 


5.3,5 Twist constraints 

The dominant design variable in most optimization computations is the twist gradient, 
TWG(iseg), over each segment, iseg. The total twist at the outboard end of iseg, 

iseg 

TWR(iseg + 1) = X TWG(jseg) * SL(jseg) (5-33) 

jseg=l 


where SL(jseg) is the segment length. Since TWR is an accumulative function of the 
TWG variables lying inboard it is possible that excessive changes in TWR occur, 
particularly as one approaches the tip, rendering the linearization assumption invalid, or 
destabilizing the subsequent wake analysis. In order to prevent such changes inequality 
constraints are adjoined to the optimization analysis limiting the twist change according 
to: 

iseg 

ATWR(iseg + 1) = X ATWG(jseg) * SL(jseg) £ 0.3° (5-34a) 

jsegeld 
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(5-34b) 


iseg 

ATWR(iseg + l) = £ ATWG(jseg) * SL(jseg) £ -0.3° 

Jsegel d 

where iseg ranges from l,...,nseg and belongs the set of segments, I<j={iseg: 

TWG(iseg)e d,}. i.e., if some of the TWG variables are held fixed, such as when only the 
outboard segments are allowed to vary in the analysis, then those fixed variables are 
skipped in the above summations. 


5.3.6 User supplied constraints 

Preliminary coding designed to incorporated user defined constraints has been 
developed (but not yet tested). It takes the form of a subroutine which allows the user to 
define inequality constraints and equality constraints as functions of the design variables, 
dj. The intended goal is to permit the incorporation of constraints derived via external 
analyses such as other forward flight studies or structural dynamics considerations. It is 
assumed that the functional dependence of these constraints upon the dj can at least be 
approximated, say by polynomial curve fitting, so that the constraints and their first 
derivatives can be defined based upon the information supplied to the subroutine. In 
principle, cost functions could also be user defined and used to augment one of the 
KOBJ(l) options. Instead, the anticipated manner of incorporating alternative cost 
functions will be limited to constructing inequality constraints such as those in 
Eqs. (5-24) and (5-25). This follows from an assumption that accurate evaluation of these 
external cost functions is better accomplished using independent models, but that their 
inclusion in a EHPIC/HERO optimization might be desirable to prevent designs yielding 
unacceptable values for these cost functions. 


5.3.7 Summary of optimization formulation 

The general optimization problem that may be formulated within EHPIC/HERO is 
summarized: 

Minimize j Cost = J + J x Ax + J y Ay + Jj Ad J 
subject to: 

Q £ gj + Qix Ax + Q iy Ay + Ad 
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The subscripts (•), and (*) c refer to inequality constraints and equality constraints (beyond 
those associated with the wake) respectively. These include constraints upon the twist, 
Eq. (5-34), and any constraints input by the user. 

If the user sets KOPT£3, then one solves for (Ajjc Ay }: 



(5-36) 


using Gaussian elimination. The optimization task is transformed by substituting 
{Aic Ay } back into the remaining equations: 

Minimize { Cost = J + J^* Ad) 


subject to: 
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(5-37) 



5.4 Numerical Considerations 

For certain options of KOPT, the simplex algorithm and BFSG update procedure can 
involve large arrays containing several hundred rows and columns. For problems of this 
size, numerical accuracy considerations become an important issue. Some of the steps 
taken in this regard are summarized in Section 5.4.1. Another consequence of such large 
problems is the increase in the computational time required to both compute the 
constraint derivatives and to perform the pivoting operations in the simplex scheme. 
Thus, various options that allow some of these computations to be skipped have been 
explored and are presented in Section 5.4.2. Finally, convergence tests are given in 
Section 5.4.3. 

5.4. 1 Treatment of round-off error 

An important consideration in implementing the SLP and SQP algo rith ms is the 
presence of round-off error. This is particularly true of larger problems involving several 
hundred constraints (e.g., incorporation of the equality constraints corresponding to the 
equations in the EHPIC/HERO code relating the wake point positions, velocities, etc.). 
Several steps were taken in order to attenuate or else accommodate the effects of round- 
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off error. Prior to the optimization analysis, the simplex tableau is preconditioned by 
rescaling all entries so as to be of comparable magnitude. First, for each row the entry 
with the maYimnm absolute value is found and the row divided by that entry. The same 
is then repeated for each columns. The latter operation is tantamount to rescaling the 
state variables, xj, and is used in determining corresponding scaling factors. Thus, 
denoting the simplex tableau entries by ajj, the scaling operations proceed as: 


a jj 

a ij 4_ T’ ’ 
*1 

rj = max{ 1 a,j 1 j 

[ , i=l,...,nrow 

(5-39a) 

ajj 

C J 

Cj = max j 1 ay 1 j 

| , j=l,...,ncol 

(5-39b) 

where for linear programming. 



nrow = m+1 

; ncol=n+l 


(5-40a) 


and for quadratic programming, 

nrow= m+n+1 ; ncol=2n+m+l (5-40b) 

The factor, 1/cj, is a representative scale for variable, xj, and is utilized elsewhere in the 
optimization algorithm when conducting comparison operations and convergence tests. 

The use of double precision for the complete tableau has the undesirable consequence 
of a four-fold increase in memory requirements. Therefore the use of double precision is 
restricted to only certain key variables and parameters used in the accumulative summing 
operations. In particular, double precision is employed for the auxiliary cost and reduced 
cost coefficients used in finding a feasible solution and, for quadratic programming 
problems, in satisfying the Kuhn-Tucker. Finally, consideration of round-off error and 
tableau dimension is made when specifying the tolerance and convergence parameters 
used in the optimization algorithm. 


5.4.2 Efficiency improvement options 

The additional computation that must be performed for optimization is categorized 
into 1) calculation of the influence coefficients pertaining to the design variables, and 2) 
execution of the SLP or SQP routine to determine the improved design. In many cases, 
the various influence coefficients change slowly with the design. Furthermore, unless the 
design is near-optimal, the design changes at each optimization step are the same, i.e., 

AX/IIAXII is constant. These fortuitous properties can be put to good use in diminishing 
the computational cost By taking larger steps in design space between optimization 
calculations and design related influence coefficient evaluation, the total number of 
operations is substantially reduced, with an accompanying trade-off of a possible 
lessening in robustness. Rather than effecting the full design change immediately after 

AX has been determined, this design change is implemented over a series of fractional 
steps. The wake analysis is reconverged during each fractional steps. This improves the 
robustness of the overall computation since the linearity assumptions implicit in the 
iterative solution in the wake analysis are more likely to be valid over one fractional step 
than over the full design change. To this end, the user supplies the input parameters, 
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KWKN and NFRAC, where NFRAC is the number of fractional steps and KWIKN has 
the following meanings: 

KWIKN=0 -No modification of the optimization procedure. The influence 
coefficients and the optimization calculation are executed every 
time the wake analysis is re-converged. 

KWIKN=1 - The influence coefficients pertaining to the design variables are 
updated once every NFRAC steps, but the optimization routine is 
executed every time the EHPIC/HERO wake analysis has 
reconverged. 

KWIKN=2 - After every NFRAC steps, the influence coefficients corresponding 
to the design variables are re-evaluated. The move-limits are 
expanded: 

Ax'"** <- NFRAC * Ax™” 

and the optimal design change, AX. computed. The design change 
is then implemented over the next NFRAC steps as: 

1) k=l 

2 > Xk «- Xk-l + nmac « 

3) Reconverge wake analysis 

4) k=k+l 

5) If k < NFRAC go to step 2 

The simplex tableau is stored in the form suggested by Kiinzi (Reference 57) which 
affords a m by m reduction in tableau size for the LP case and (m+n) by (m+n) reduction 
in the QP case. An by n block can be eliminated in the QP situation by recognizing that 

the tableau coefficients for the Lagrange multiplier vectors, jii and j!2> a re identical 
except for sign. This is due to the twin observations that: 1) the initial tableau entries 

corresponding to pi and (J.2 are simply [-I n ] and Pn] respectively so that the above 
relationship holds at the start of the simplex manipulations; and 2) the pivoting operations 
basically consist of row manipulations where each entry in a row is operated upon in the 
same manner. Thus rows may be multiplied by the same scalar, and scalar multiples of a 
row may be added to another. It thus follows that columns that arc multiples of each 
other retain that property subsequent to the pivot operations and that furthermore, the 
constant of proportionality remains unchanged. 


5.4.3 Storage requirements 

For SQP optimization, storage space on the order of (2n+m) by (n+m) is needed 
where n is the total no. of variables under consideration (including in addition to the 
vector of design variables, the vector of bound circulation perturbations and normal/bi- 
normal wake point perturbations unless these are explicitly eliminated via Eq. (5-36); and 
m is the total no. of equality and inequality constraints imposed upon the problem (not 
including upper and lower bounds imposed upon the variables). The required storage 
arises from the fact that there is a total of 3n+m variables comprised of the n state 

variables, xj, 2n Lagrange multipliers, Hu and pji. associated with each of the lower and 

upper bounds of xj, and m Lagrange multipliers, Xj, associated with the imposed 
constraints, gj; and m+n constraints including the m original imposed constraints and the 

additional n Kuhn-Tucker constraints. The simplex tableau columns corresponding to (i^ 


78 



and H 2 i differ only in sign and so only one of the columns must be stored (see Appendix 
A). This reduces the tableau column size from (3n+m) to (2n+m). 

The SLP approach requires storage space on the order of n by m. Both the SLP and 
SQP schemes implicitly enforce the upper and lower bounds on xj as explained in 
Appendix A. Thus elimination of the 2n constraints required during the Phase I effort to 
enforce these bounds are no longer necessary. 


5.4.4 Convergence criteria 

The selection of suitable convergence tolerances has proved to be a non-trivial task in 
the multi-dimensional optimizations considered here. The different scales for the design 

variables (and the Xc and 3 ) are determined in the EHPIC/HERO code and taken into 
account in convergence tests. However, the chief obstacle in deriving generally 
applicable convergence tests is due to the large disparities in gradients and curvatures for 
different cost functions and designs. For example, the thrust is typically one or two 
orders of magnitude more sensitive to changes in the twist gradient than to sweep 
perturbations. Thus it can happen that the twist distribution has essentially converged to 
the optimal design, while the sweep distribution is still undergoing significant 
modification. A test based upon changes in cost per step might cause premature 
termination of the optimization since the changes in thnist level are 'small' in comparison 
to the changes that occurred while the twist distribution was still changing significantly. 
Thus consideration of the local gradients should be factored into the convergence tests. 
Other elements that complicate the definition of an appropriate convergence criterion 
include the possibility of inaccurate approximations to the Hessian matrix, [B], in SQP 
(can be countered by resetting [B] near a prospective minimum), possible discontinuities 
in some of the first order derivatives which can arise in implementing a stall model, and, 
importantly, round-off error. 


The following 'classical' convergence tests are available in EHPIC/HERO and are 
based upon the user supplied tolerance parameter, TOL. The optimization terminates if 
any of these tests are satisfied. If the user deems the optimization to have terminated 
prematurely, then TOL may be reduced and the restart option of EHPIC/HERO invoked. 


a) 

b) 

c) 


I JQ - JQ-l I < TOL * Min. { I J° I , R ) , for 3 successive optimization steps 
I jq - jq-1 1 < TOL * Min. { I JQ I , R } , for 3 successive optimization steps 


max.-j 



SC(i) 


} £ TOL 


d) 

c) 



Axj 113 * 

SC(i) 


q > qmax 


< TOL ( for SLP optimization only ) 
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Here, J4 is the objective (or cost) function at optimization step q, SC(i) is a representative 
scale computed in EHPIC/HERO for the design variable under consideration, x- 1 is the 

value of the i-th design variable at optimization step, q. Ax}" 8 * is the move limit for the i- 
th design variable and is initialized to DSDMAX at q=0, qmax is the m a xim u m allowable 
number of iterations, and the parameter, R, is an estimate of the reduction in cost that can 
be achieved for the given move limits, DSDMAX(i), and is computed from 


dJ 


R = 0.2 * Y * DSDMAX(i) 

i i 


(5-41) 


As explained in section 5.2.4, when employing SLP the move limits, Ax-" 2 *, are reduced 
when proximity to the optimum is detected. This process continues until criterion d) is 
satisfied. 
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6. RESULTS OF SAMPLE PROBLEMS: ROTOR PERFORMANCE 


6.1 Previous Validation Work 

As noted in the discussion above, meaningful design optimization requires that 
the methods used to model rotor performance produce accurate results. Extensive 
performance correlation studies were carried out during the development of the original 
EHPIC Mod 0.0 code (Ref. 9), and additional studies were executed by NASA personnel 
(Ref. 6) as well as by users of the EHPIC (Mod 1.0) version (Ref. 7). These correlation 
studies involved tests of many different rotor configurations, including two-, three-, four-, 
five- and six-bladed helicopter main rotors, three tiltrotor configurations, and several tail 
rotors. A wide range of designs were examined including tapered planforms, swept 
planforms, twisted planforms and combinations of each. Nonlinear as well as linear twist 
schedules were investigated including the very high twist levels characteristic of tiltrotor 
blades. Several rotors in ground effect have also been studied (Ref. 10, 11). Finally, 
recent studies (Ref. 60) of more limited scope have shown that EHPIC compares quite 
favorably and in some res pects improves upon the performance results obtained using 
much more CPU-intensive CFD analyses. 

These studies have concluded that the EHPIC code produces very good 
performance prediction across a broad range of rotor designs. As with any numerical 
analysis, the results obtained are sensitive to certain key input parameters. However, the 
results described in Reference 7 were particularly encouraging on this point, in that 
calculations carried out over several very different rotor systems produced consistently 
good correlation for a fixed set of inputs - i.e., no use of numerical "dials" was necessary 
to obtain good results in particular cases. One of the objectives of the present study was 
to ensure that the changes implemented during the development of the HERO analysis 
produced the same consistently good performance correlation as its parent EHPIC code. 
Some changes in the results were unavoidable given the necessity of restructuring the 
code to adapt to the requirements of carrying out performance optimization. Also, as 
previously discussed, several improvements in the aerodynamic model have been 
implemented to eliminate particular limitations of the baseline analysis. The correlation 
studies that follow are intended to illustrate the fundamental consistency of the 
optimization hover performance model with the original EHPIC code as well as to 
highlight the improvements made. 

In each of the calculations below, the EHPIC and EHPIC/HERO performance 
computations used the same number of trailing filaments, extent of free and prescribed 
wake, vortex lattice discretization, and core modeling. To the extent that the results of 
the two analyses differ, the differences are analyzed in each particular case. Unless 
otherwise specified, standard atmospheric conditions are assumed corresponding to a 
speed of sound of 340 m/s (1117 fps) and atmospheric density of 1.205 kg/m 3 (.002378 
slugs/ft 3 ). Since NASA is currently operating the Mod 1.0 variant of EHPIC, 
comparisons of HERO and EHPIC here are made with the 1.0 version. 

6.1.1 NAS A/NACA test rotors 

A data set from an NACA rotor tested in the 1950's provides a simple 
configuration with which to start the present correlation study. The first test rotor 
examined was that described in Reference 61, which described integrated performance 
results on a two-bladed rotor with a radius of 8.17 m. (26.8 ft.) and a constant chord of 
0.58 m. (1.91 ft.). The blades featured 8 degrees of linear twist and were tested at tip 
Mach numbers of 0.28 and 0.66. An NACA 0012 airfoil section was used across the full 
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span of the blade. The wake model here used six filaments trailing from the rotor blade. 
Since the test was conducted on a whirl stand elevated 1.5R above the ground, the ground 
plane model was invoked to correct for the effect of the image system of vortices. The 
blade model used 30 vortex quadrilaterals along the span and one chordwise. 


Figures 6-1 and 6-2 show the correlation achieved for both tip Mach number 
cases. The correlation is good across most of the range tested, though some differences 
start to appear at high thrust. Similar results were obtained with the original EHPIC 
code, as described in Reference 9. It is instructive to compare the details of the results 
obtained with EHPIC Mod 1.0 against those obtained with EHPIC/HERO since some 
modifications have been made to the basic performance model. For simplicity, a case 
with Map = 0.28 is chosen. Figure 6-3 shows the results of the integrated power and 
thrust results for runs using EHPIC Mod 1.0 and EHPIC/HERO over a range of 
collectives. The predictions of the two codes are very close, but some small differences 
do appear to be present 


Further details are provided by a close examination of a single case. For a root 
collective pitch angle of 14 degrees , the following integrated performance results were 
obtained: 


EHPIC/HERO 


EHPIC Mod 1.0 


C T 0.003239 

Cni 0.0001449 

Cop 0.00004454 

Cq 0.0001895 


0.003261 

0.0001429 

0.0004446 

0.0001875 


These results indicate the typical size of differences in performance predicted by 
EHPIC/HERO and EHPIC Mod 1.0. The profile power calculation yields essentially 
identical results, as would be expected given that the same airfoil drag characteristics 
were read in, though in a different format in the new code. Differences in the thrust and 
induced power are larger, and are attributable primarily to shifting the vortex filament 
release points to the trailing edge as described in Section 2.3.1. This change was 
necessary for the implementation of the high-resolution wake roll-up calculation 
described in Section 3. Figure 6-4 shows the distributed thrust, induced power, and total 
power for the two calculations. Most of the differences between the two calculations 
appear in the tip region, while the computed loading inboard of roughly 0.85R is nearly 
identical. 


6.1.2 CH-47B rotor 

Performance tests for the CH-47B main rotor are described in Reference 62. The 
rotor has three blades, each with a constant chord of 0.64 m. (2.1 ft) and a radius of 9.14 
m. (30.0 ft.). The blades feature -9.14 degrees of linear twist and use a 23010-1.58 airfoil 
section. The data in Reference 62 was taken at rotor rotation rates of 230, 240, and 288 
rpm. The case examined here is 240 rpm, corresponding to a tip speed of 229.6 m/sec 
(753 fps). 

The wake model used six filaments trailing from the span, with roughly 2.5 turns 
of free wake followed by an additional 1.5 turns of prescribed wake. Thirty vortex 
quadrilaterals were used across the span with one chordwise. The integrated performance 
predicted for this case is shown in Figure 6-5, indicating good agreement up to high 
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Figure 6-1. Predicted and measured performance of the two bladed NASA 
TN4357 rotor at Mtip = 0.28 . 
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Figure 6-2. Predicted and measured performance of the two bladed NASA 
TN4357 rotor at Mtip = 0*66 . 
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Figure 6-3. Comparison of predicted performance using the EHPIC Mod 1.0 and 
EHPIC/HERO codes for the NASA TN4357 rotor. 
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Figure 6-4. Predicted distributed loading for the NASA TN4357 rotor. 
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Figure 6-5. Predicted and measured performance for the CH-47B main rotor. 
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Figure 6-6. Predicted and measured performance for the UH-60A model rotor. 
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thrust levels with some deviation beyond that point. The performance predictions are 
similar to those obtained with EHPIC Mod 0.0 in Reference 9. 

6.1.3 UH-60A rotor 

Recently acquired data on an extensively instrumented UH-60A model rotor has 
provided a wealth of potentially useful correlation information for study using 
EHPIC/HERO. References 63-65 that describe this test were acquired only shortly 
before the end of this effort, and analysis and interpretation of the results is still ongoing. 
However, certain major results of the computational studies done to date are available and 
are presented here. 

The model tests described in References 63-65 involved a 8.86 m. (9.4 ft.) 
diameter, Mach scaled model of the four-bladed UH-60A main rotor. The model featured 
a chord of .092 m. (0.304 ft.), and was operated at a range of tip Mach numbers between 
0.55 and 0.70. This planform features 20 degrees of sweep starting at the 92% radial 
station, as well as a distinctive and highly nonlinear twist distribution that includes a twist 
l>ucket' over roughly the outer 10% of the blade span (the UH-60 twist distribution is 
discussed further in Section 7). The particular case considered here is at tip Mach 
number 0.65. The blade uses an SC1095 airfoil across the entire span, except for the 
'droop nose' SC1095R8 airfoil between 0.47R and 0.82R. The computational model uses 
45 equally spaced quadrilaterals spanwise and one chordwise. Six free filaments are 
used, with one capturing the wake of the tip region and five used to model the inboard 
sheet 


Figure 6-6 shows the prediction of integrated loading for Run 68 of Reference 63. 
The power prediction is close over most of the range surveyed, though with some 
underprediction of the power at higher thrusts. The prediction of thrust as a function of 
collective pitch is likewise close, as shown in Figure 6-7. At high thrust levels, 
EHPIC/HERO tends to slightly overpredict the actual thrust 

The tests described in References 63-65 involve some runs at very high thrust 
coefficient, so it was judged a suitable test case for the static stall model discussed in 
Section 2. However, the range of thrust coefficients tested did not in fact pass the cin, M 
level for the airfoil tables tested, as indicated in Figure 6-8. The computations do show 
that the qualitative behavior of the present model is reasonable, predicting a drop-off in 
thrust growth with collective (Figure 6-8) along with a rapid increase in power required 
(Figures 6-9 and 6-10) once the stall model is activated. (The primary purpose of the stall 
model is to limit performance at twist angles above the section stall angle during the 
optimization process. Most current planform designs operate far below the stall regime 
so the stall model will not affect comparisons with test results.) 

As noted above, References 63-65 contain a wide range of experimental 
measurements suitable for further correlation studies, though the recent (early 1992) 
release of the data has precluded detailed consideration here. In particular. Reference 65 
documents measured structural deflections that could be of considerable value in 
validating the structural deflection model described in Section 4. It is anticipated that 
future validation efforts will include such correlations, as well as studies of the spanwise 
loading and wake geometry data contained therein. 

6.1.4 XV-15 rotor 

As noted in the wake geometry calculations in Section 3, the XV-15 rotor is a 
three bladed configuration with a radius of 3.81 m. (12.5 ft.) and a constant chord of 
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Figure 6-8. Rotor thrust predictions versus collective with and without stall 
model (UH-60A model rotor). 
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Figure 6-9. Rotor profile torque predictions versus collective with and 
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0.354 m. (1.16 ft.). The twist distribution features a nonlinear twist distribution with a 
total washout of approximately forty degrees across the span. Each blade uses five 
NACA 64-class airfoils across the span. The tip Mach number for the cases studied here 
was 0.69, and the measured performance is drawn from Reference 45. The same blade 
model with 45 vortex quadrilaterals span wise and 1 chordwise used in Section 3 was used 
here, though the present calculations employ eight equally spaced filaments to obtain 
good resolution of the spanwise distribution of induced velocity. 

Figure 6-11 shows the correlation achieved across the moderate to high thrust 
range of most practical interest The calculations are compared to two different runs from 
Reference 45 that bracket most of the range of data taken in the test. The correlation is 
good across the range examined, though there is some tendency to underpredict the power 
at low thrust levels. 

6.1.5 V-22 rotor 

Reference 27 discusses performance measurements taken on a 0.658-scale V-22 
rotor. Like the XV- 15, the V-22 is a three-bladed rotor with a large, nonlinear washout 
across the span, with a total of approximately 40 degrees of twist across the span. The 
scaled V-22 rotor has the same radius as the XV-15 but also features a tapered planform, 
with the chord ramping down from 0.58 m. (1.9 ft) at the root to 0.405 m. (1.33 ft) at the 
tip. The tip Mach number for these tests was 0.68, corresponding to a tip speed of 231.7 
m/sec (759.7 ft/sec) at standard conditions. The wake and blade model were identical to 
those used for the XV-15 case. The results obtained here were likewise similar. Figure 6- 
12 shows the correlation of integrated performance, again indicating good general 
agreement, though with some errors appearing at both the low and high end of the range 
examined. 

6.2 Performance in Axial Flight 

Only limited measured performance data is presently available for proprotors in 
high speed axial flight. One of the designs of high current interest in this area is the V-22 
rotor, whose performance is discussed in Reference 66. Felker describes the following 
closed-form semi-empirical expression for integrated rotor power that has been quite 
successful in capturing on-design propulsive efficiency for the V-22 (Ref. 67). 


Cq = Cqp + pC T + ^ 


CQp = ^C^Qo 


(l+^ 2 )Vl+li 2 +^t 4 ln 


h +VT 




where f c is an empirically-derived compressibility adjustment factor for the profile drag: 

\ 


fc = 


1, if M .75 £ 0.63, otherwise 
\l + 42.51 (M.75 - 0.63 ) + 3476 ( M.75 - 0.63 J 4 f 

M ” - ^ 4 ^ 
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Figure 6-11. Predicted and measured performance for the XV- 1 5 rotor. 
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In order to test the capability of EHPIC/HERO to successfully analyze rotors in high 
speed axial flight this formula was used in place of directly measured performance. 

Curves of rotor propulsive efficiency pCf/Cp as a function of thrust coefficient were 

made for two flight conditions: first, p = 0.34 and Mhd = 0.62 and, second, p = 0.67 and 

Mjiei =0.59. 

Since rotor wake effects are very weak in high speed axial flight, very simple 
wake models are appropriate for these flight conditions. In this case, five free filaments 
with one turn of free wake were used, though the filaments are quickly convected 
downstream, behaving essentially as a kinematic wake. With thirty vortex quadrilaterals 
across the span, good resolution of spanwise loading could be expected. The integrated 
performance results shown in Figures 6-13 and 6-14 indicate that the major features of 
rotor performance are in fact being captured. For the moderate speed case shown in 
Figure 6-13, the prediction accuracy is good across most of the range, while for the high 
speed case shown in Figure 6-14, a consistent underprediction of propulsive efficiency is 
apparent. The presence of such a constant decrement suggests inadequacies in the profile 
power calculation, possibly in the high Mach number section data for the airfoils used on 
the rotor. 

Figure 6-15 shows a comparison between predicted and measured performance 
for the V-22 rotor in low speed axial flight A small range of EHPIC/HERO results are 

shown at three advance ratios, p = .0236, .0313, and .0381. Only one test data point was 
available at each of these advance ratios. EHPIC/HERO shows good agreement at the 

two lower values but begins to over-predict the torque at p = .0381. Further test data will 
have to be obtained in order to determine whether this trend occurs with increasing 
advance ratio or increasing thrust 
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Figure 6-13. EHPIC/HERO prediction of V-22 propulsive efficiency in high speed 
axial flight: advance ratio 0.34, Mhel = 0-62. 
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Figure 6-14. EHPIC/HERO prediction of V-22 propulsive efficiency in high speed 
axial flight: advance ratio 0.67, Mh c i = 0.59. 
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7. RESULTS OF SAMPLE PROBLEMS: DESIGN OPTIMIZATION 


This section presents results that illustrate the application of the EHPIC/HERO code 
in optimization calculations for rotors in hover and axial flight. The initial focus of 
interest is on the problem of power minimization at a given rotor thrust for rotors in 
hover. This is the rotary-wing analog to the classical fixed-wing drag minimization 
problem that still commands attention today. The power minimization problem is 
likewise critical in rotorcraft design because of its importance in sizing the powerplant to 
satisfy payload requirements. Here, we analyze several representative configurations to 
both illustrate the capabilities of the present formulation and to gain insight into 
promising general design strategies. 

As noted in Section 5, EHPIC/HERO has the ability to define a variety of objective 
functions besides total rotor power. Gross thrust, hover Figure of Merit, axial propulsive 
efficiency, and individual components of total power (i.e., induced and profile) may be 
selected, as may weighted combinations of these quantities. We will show results of 
selected problems executed using some of these alternative objective functions for 
hovering rotors, including thrust maximization for a fixed power. 

EHPIC/HERO is also suited to the analysis of rotors in axial flight, as discussed in 
Section 6. The design requirements of rotors in high speed axial flight (e.g., tiltrotors in 
cruise) are distinct from those of rotors in hover. The influence of the wake is much 
reduced, so much so that uniform inflow approximations can often be used with some 
success, as the discussion of performance results in Section 6 indicates. However, even 
though the free wake feature of the present analysis is not as important as in hover, the 
lifting surface/vortex lattice aerodynamic model in EHPIC/HERO can be used to good 
effect here, since the design evolution of proprotors often calls for swept and variable 
chord planforms. The performance results for axial flight calculations to date have been 
promising, and it is anticipated that application of the design optimization process to 
rotors in high speed axial flight will yield informative results. 

Finally, it should be noted that the purpose of these demonstration calculations was 
not to comprehensively exercise the extensive set of design optimization options within 
EHPIC/HERO. The design evolution in particular cases is a function of the constraints 
imposed, the algorithms selected, and to some extent of the numerical resolution of the 
blade and the wake that is compatible with the user's computational constraints. The 
immediate objective here is to present calculations that illustrate some of the major 
capabilities of the present analysis. As will be discussed, certain broad trends were 
observed in the calculations performed to date that appear to represent generally desirable 
design strategies for rotor blades, while other cases provoke as many questions as they 
answer. 

7.1 Sample Calculations in Hover: Rotor Power Minimization at Constant Thrust 

7.1.1 Conventional low- twist helicopter designs 

Demonstration calculations of performance optimization were carried out on a 
planform similar to the UH-60A model rotor analyzed in Section 6, but without the tip 
sweep and twist 'bucket' characteristic of that blade (Figure 7-1). The resulting blade was 
straight and untapered, with a radius of 4.68 ft and a constant chord of 0.303 ft. The 
resulting twist distribution yields a -16 deg. twist rate across the blade radius. An initial 
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Figure 7- 1 . Twist distribution of the UH-60A model rotor blade. 
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Figure 7-2. Figure of merit history during design optimization of a four-bladed 
UH-60A-class rotor: twist and tip sweep. 
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Figure 7-3. Initial (solid) and final (dashed) bound circulation distributions for 
the twist/tip sweep optimization. 
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condition with thrust coefficient of 0.00664 was selected, which corresponds to Cj/ct of 
0.0805 and a Figure of Merit of 0.721. 

The same wake configuration used in the performance studies of the UH-60A 
model blade in Section 6 was used here. The blade was discretized into ten unequal 
segments, ranging in width from 0.1 25R to 0.025R, with the smaller elements 
concentrated near the tip to capture the expected rapid design gradients in that region. 
Unless otherwise noted in the text below, the sequential linear programming (SLP) option 
was used, and all blades were considered rigid (no structural deformation allowed). As 
described previously, while the structural model has been validated through the use of 
idealized calculations and comparisons with analytical results, a calibration against 
measured blade deflection has not yet been performed. 

Several different combinations of design variables were selected for investigation. 
These were: twist and tip sweep; twist and tip anhedral; and twist and variable chord. In 
each case, the twist was constrained to have twist changes between -8.5 and +1.5 degrees 
on each segment. This 10 degree range was judged to give ample latitude for potentially 
interesting design trends to appear without producing unrealistic configurations. 

Twist and tip sweep : In this case, the outer 7% of the blade was allowed to sweep 
up to +/-30 deg. The results of this case are summarized in Figures 7-2 to 7-4. Figure 7- 
2 shows the history of the evolution of the hover Figure of Merit during the calculation, 
as it rises to roughly 0.75 from its initial value of 0.721. Though the calculation dithers 
considerably around the plateau, the bulk of the improvement comes in the first fifteen 
design optimization steps, during which time the tip sweep angle ramps up monotonically 
from zero to the upper bound of 30 deg. The twist distribution undergoes dramatic 
evolution during this period also, as shown in Figure 7-4. The distribution begins to 
develop the twist bucket near the tip originally built in to the UH-60A. Some tendency to 
increased twist over the inboard part of the blade is also observed. The overall tendency 
of the evolution of the bound circulation distribution, though, is to flatten out the peak in 
circulation seen near the tip and to produce a much more uniform distribution (Fig. 7-3). 
As will be seen, this trend is evident in a wide variety of die cases studied to date. 

Twist and tip anhedral : This case produces results in many respects qualitatively 
similar to the previous exercise. The same discretization of the blade was used, with ten 
segments of decreasing span being used. Here the outer 17% of the blade was allowed to 
deflect in anhedral, with the deflection limited to +/-15 deg. Again the twist distribution 
was allowed to vary across the full span. Figure 7-5 shows the history of the Figure of 
Merit as a function of number of optimization steps. Once again, a rapid rise is observed 
during which the anhedral angle (tip droop) ramps rapidly to its prescribed limit. The 
relatively large fluctuations observed in the Figure of Merit in steps 10-15 are due to 
adjustments in the strength and position of the wake. These changes settle down to a 
relatively stable plateau after steps 15-20. The predicted twist distribution associated 
with the drooped tip is shown in Figure 7-7, showing a smaller twist bucket near the tip, 
though again with larger twist angles toward the root The resulting bound circulation 
distribution at the termination of the calculation contains an anomalous trough that 
appeared to be diminishing as the calculation plateaued around step 40 (Figure 7-6). 

Twist and chord: Here, both the chord and the twist distribution were allowed to 
vary. The chord was constrained to stay at its baseline value of 0.305 ft over the inner 
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Initial (solid) and final (dashed) geometric twist distributions for 
the twist/tip sweep optimization. 
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Figure of merit history during design optimization of a four-bladed 
UH-60A-class rotor: twist and tip droop (anhedral). 
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Initial (solid) and final (dashed) bound circulation distributions for 
the twist/tip anhedral optimization. 
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Figure 7-7. Initial (solid) and final (dashed) geometric twist distributions for 
the twist/tip anhedral optimization. 
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Figure 7-8. Initial (solid) and final (dashed) bound circulation distributions for 
the twist/chord optimization. 
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40% of the span, while outboard of this point increases of up to 33% over the baseline 
were allowed; the minimum permissible chord in this region was set at 0.1 ft, or one 
third of the baseline. In this case, the Figure of Merit increased from the baseline value 
to just under 0.77 in thirty steps after which point the performance leveled off. The chord 
distribution tapered gradually down to 0.18 ft. from the end of the fixed-chord segment, 
while chords as large or larger than the baseline appeared in the immediate vicinity of the 
tip. Figure 7-9 shows the twist distribution obtained for this case, which is similar to the 
other two cases just discussed. The tendency to a constant circulation distribution is 
shared with the previous calculations. Further investigation will be required to determine 
if the increased chord at the planform tip is a consequence of the lattice discretization 
used or if it represents a physically meaningful event 

These calculations indicate that representative baseline configurations can achieve 
substantial performance increments of 3-4 points in Figure of Merit without excessive or 
radical design departures. What is of more interest than the particular numerical results 
achieved is the appearance of two trends in the design evolution: the tendency to 
smoother, flatter bound circulation distributions and the appearance of features 
resembling the twist bucket shown in Figure 7-1 near the tip. Some indications of the 
desirability of flatter circulation distributions were present in the sample cases studied in 
Reference 4, and the trend is reasonable in light of classical studies of the advantages of 
uniform downwash fields. A flat circulation distribution minimizes the strength of the 
wake circulation and hence, in general, minimizes the induced power caused by wake 
downwash on the blade. 

The appearance of the twist bucket, though found empirically desirable by the 
designers of the UH-60, is not a feature that emerges naturally from simple performance 
analysis. It is noteworthy that the trend toward this particular design feature is driven 
largely by the minimization of induced rather than profile power. Repeating the 
calculations shown above with the profile power artificially set to zero - while obviously 
leading to lower overall power - still yields twist distributions with the down-up 
distribution near the tip characteristic of the bucket in Figure 7-1. A likely interpretation 
is that the down-up twist distribution contributes to the leveling of the circulation 
distribution by dropping the load near the peak and increasing it at the blade tip. A 
similar phenomenon can be seen near the root 

7.1.2 Tiltrotor 

The design optimization of a tiltrotor in hover presents a substantially different 
challenge from the computations just described for more conventional low-twist baseline 
designs. A representative case was examined to investigate possible design improvements 
in baseline tiltrotor designs. The case considered here for illustration is an XV- 15 rotor 
with the same operating state and planform as those studied in Section 6. For this case, 
optimization of twist and tip sweep was considered. The twist changes were effectively 
unconstrained, while sweep was constrained to be zero except for the last 10% of the 
blade, which could adopt a sweep of +/-30 deg. The thrust coefficient of the base case 
was 0.0127, while the Figure of Merit was 0.789. 

The evolution of the Figure of Merit is shown in Figure 7-10. The history is 
notably uneven, but levels off around 0.815, corresponding to roughly a 3% power 
reduction at constant thrust The unevenness in the advance of the Figure of Merit is 
attributable in part to oscillations in the tip sweep angle as it ramps up from zero to 30 
degrees at step 30. The total performance increment for this case is relatively modest but 
this is not surprising if the initial and final bound circulation distributions are considered 
(Figure 7-11). The initial bound circulation distribution here was already relatively 
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Figure 7-10. Figure of merit during design optimization of an XV- 15 tiltrotor. 
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Figure 7-11. Initial (solid) and final (dashed) bound circulation distributions for 
the XV- 15 power minimization. 



Figure 7-12. Initial (solid) and final (dashed) geometric twist distributions for the 
XV- 15 power minimization. 
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uniform and left little latitude for subsequent modifications. It is noteworthy, though, 
that the principal change made by the analysis is the introduction of an up-twist near die 
blade tip (Figure 7-12). This feature is qualitatively similar to that found in the 
calculations in Section 7.1. Its appearance in a configuration radically different from the 
low-twist, four-bladed design just discussed suggests that the introduction of this type of 
twist distribution may be a genetically desirable design feature for a wide variety of 
rotors. 

7.2 Alternate Objective Functions 

All of the cases considered to this point have been power minimizations at constant 
thrust. Many other types of performance optimization problems are of practical interest, 
and a wide range of objective functions can be accommodated, as discussed in Section 5. 
A sample calculation was undertaken to demonstrate this capability; the problem involves 
thrust maximization at constant power, corresponding conceptually to the task of 
maximizing payload for a given power plant. The base case used here was the same UH- 
60A-class rotor studied in Section 7.1. 

Starting from the baseline thrust coefficient of 0.00664, the design was allowed to 
evolve with the same constraints imposed in the twist and tip sweep optimization 
described in Section 7.1. The resulting histories of thrust and Figure of Merit are shown 
in Figures 7-13 and 7-14. A thrust increase of roughly 2.9% is achieved, with a 
corresponding increase in Figure of Merit of approximately .038 as a consequence of the 
constant power constraint. The spanwise loading and geometric twist (Figures 7-15 and 
7-16) for this 'thrust optimized' design are similar to the distributions obtained in the 
'power optimized' design discussed above. 

As described in Section 5, many other candidate objective functions can be 
accommodated within EHPIC/HERO, including weighted functions of thrust, power (and 

its induced and profile components), Figure of Merit, and propulsive efficiency (jiCj/Cp) 
for axial flight The latter case is now considered. 

7.3 Axial Flight: Tiltrotor/Proprotor Case 

Existing and proposed tiltrotor designs call for blade designs that can operate 
efficiently at very high axial flow rates, typically as high as 450 fps. The discussion in 
Section 6 included predictions of the performance of a V-22 tiltrotor in two cruise 
conditions. Here, a study was undertaken to investigate the design trends when improved 
designs were sought for a representative high speed case. 

Adopting a strategy analogous to that described for conventional rotors in Section 
7.1.1, the baseline configuration selected was a nonspecific but representative tiltrotor 
planform, having characteristics similar to the 0.658-scale V-22 rotor and the XV- 15 
rotor studied in Section 6. The design featured a three-bladed 3.81 m. (12.5 ft.) radius 
rotor with a constant chord of 0.457 m. (1.5 ft.) yielding a solidity of 0.088 and -40 
degrees of washout, assumed to be linearly distributed from the root to the tip. Given the 
selected operating condition of 1 12.5 m/sec (369 ft/sec or 219 kts), the advance ratio was 
0.67 and the helical tip Mach number was 0.59. Each blade uses five NACA 64-class 
airfoils across the span as in the XV-15 case discussed previously in Section 6.1.4. 

The blade was initialized with 30 constant-width vortex quadrilaterals across the 
span and one chordwise. Because of the dominance of free stream convection, a 
relatively coarse wake model with a single turn of free wake on five free filaments is 
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Figure 7-13. Figure of merit history during design optimization of a UH-60A- 
class four bladed rotor: thrust maximization using twist and sweep. 
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Figure 7-14. Thrust coefficient history during design optimization of a UH-60A- 
class four bladed rotor: thrust maximization using twist and sweep. 
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Figure 7-16. Initial (solid) and final (dashed) geometric twist distributions for the 
thrust maximization calculation. 
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adequate. The span was discretized into ten segments outboard of the cutout, each 
spanning 0.09R and having -4.0 degrees of washout. 

The calculation was set up to maximize the propulsive efficiency at a constant 
thrust coefficient of 0.0085, representing a relatively heavily loaded case. All of the 
major design variables were allowed to participate in the optimization, including twist, 
chord, sweep, and anhedral. The twist was again constrained to have a maximum delta of 
between -8.5 and +1.5 across the segment. Chord was constrained not to decrease over 
the first 0.37R of the blade, but was allowed to taper down to 50% of the baseline chord 
outboard of this. Sweep and anhedral angles were kept at 0 degrees inboard of 0.46R and 
limited to +/-30 degrees outboard of this. The blades were assumed rigid, with no 
structural deformation included. 

The improvement in propulsive efficiency over 70 optimization steps is shown in 
Figure 7-17, indicating that the propulsive efficiency has leveled off at roughly 0.933, 
representing an increment of 0.023 over the baseline value. The bound circulation 
distribution shown in Figure 7-18 indicates that the familiar trend to a more uniform load 
distribution holds here as well. The geometric twist results (Figure 7-19) show relatively 
modest changes, though the up-twist near the tip characteristic of the hover results does 
not appear. 

The baseline and modified planforms are shown in Figures 7-20 and 7-21 
respectively. Sweep, anhedral and chord distributions are shown in Figures 7-22 through 
7-24. The axial flight case exhibits trends similar to those observed in the previous cases. 
The sweep and chord distributions adjust in a manner that flattens out the circulation 
distribution as much as possible which reduces the induced torque. This is best seen by 
the shape of the optimized chord distribution which is actually an inverted image of the 
circulation distribution. Though the sweep back will afford a small reduction in profile 
torque due to compressibility effects, the primary effect for rotorcraft appears to be the 
reduction of induced torque caused by reducing the circulation at the high speed tip. The 
outboard anhedral also reduces induced torque, in part by pushing the tip vortex 
downward away from the blade. Even though this latter effect is lessened in axial flight, 
the optimization algorithm will still droop the tip if any increment in performance is to be 
gained. This fact is important to take into consideration when analyzing optimized 
planforms; large excursions in sweep and anhedral often have a smaller effect on 
performance than minor adjustments in the twist distribution. It is these design variables 
having the least influence that will often change the most 

7.4 Computation Time 

Design optimization calculations are inherently computationally intensive since 
they inevitably involve repeated calls to the performance evaluation routines. As noted 
earlier in this report, the formulation of the original EHPIC code helps to reduce this 
burden since many of the influence coefficients needed to fill the tableau used to solve 
the optimization problem are computed as a matter of course in the performance 
evaluation. Nevertheless, the computational demands of EHPIC/HERO can become 
substantial as the number of design degrees of freedom are increased. The overall 
objective of the present effort in this respect was to ensure that EHPIC/HERO was no 
more CPU-intensive on 1992 computational hardware than was EHPIC Mod 0.0 on 1987 
hardware. Since EHPIC has gained acceptance in the rotorcraft industry for routine 
aerodynamic calculations, this was judged to be a reasonable criterion by which to gauge 
the usefulness of EHPIC/HERO. 
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Figure 7-17. Evolution of propulsive efficiency during the design optimization 
of a tiltrotor in high speed axial flight. 
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Figure 7-18. Initial (solid) and final (dashed) bound circulation for the tiltrotor 
propulsive efficiency calculation. 
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Figure 7-19. Initial (solid) and final (dashed) geometric twist distributions for the 
tiltrotor propulsive efficiency calculation. 
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Figure 7-20. Top view of rotor configurations for the tiltrotor propulsive 
efficiency maximization. 
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Figure 7-21. Oblique view of rotor configurations for the tiltrotor propulsive 
efficiency maximization. 
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Figure 7-22. Sweep distribution optimization of a tiltrotor in high speed axial flight. 
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Figure 7-23. Anhedral distribution optimization of a tiltrotor in high speed axial flight. 
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Figure 7-24. Chord distribution optimization of a tiltrotor in high speed axial flight. 
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This objective has been achieved. The sample calculations presented in this section 
typically run 6-10 CPU hours on an Iris 340 workstation. These run times are 
comparable to those required by complex EHPIC Mod 0.0 runs on the Micro VAX II that 
was originally used for EHPIC development during 1986 and 1987. Indeed, the present 
calculations were relatively conservative, in that they did not take full advantage of the 
options available to bypass updates of the optimization tableau. Each of the calculations 
shown in this section were run with KOPT-1, NFRAC=2, and KWIKN=2 (see Section 5 
or Reference 12 for definitions of these parameters), and it is quite likely that additional 
speed-up could have been achieved with larger steps between updates. Even so, 
additional effort to improve computational efficiency would nonetheless be desirable. 
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8. SUMMARY 


The primary purpose of this report has been to document the development and 
testing of the EHPIC/HERO rotor performance optimization code. EHPIC/HERO has 
been designed to retain the strengths of the original EHPIC free wake hover performance 
code while both extending its capabilities and adding wholly new features to compute 
structural deflection, high-resolution rotor wake flows, and, most importantly, optimized 
rotor designs for improved performance. 

The major refinements to the basic hover model include the addition of a lift stall 
model, the inclusion of scan planes for induced velocity calculations, the expansion of 
previous table look-up routines to support structural deflection calculations, and the 
revision of blade/wake coupling to facilitate both optimization and high-resolution rotor 
wake computations. The new stall model produces the correct qualitative trends for 
highly loaded rotors. The results of sample calculations of time-averaged induced 
velocities using the new scan plane capability have likewise been promising. 
Implementation of a high-resolution tip wake calculation option has also been successful, 
and testing on representative rotor configurations has demonstrated the code's ability to 
capture the structure of the rolling-up wake from blades with shallow bound circulation 
gradients near the tip. Finally, the revised approach to coupling the blade and the wake 
described in Section 2 has proved to be robust, as was found in the long performance 
optimization runs - typically involving large changes in bound circulation - described in 
Section 6. 

The formulation and implementation of a finite element model of the blade 
structure has been described. The present model uses 14-d.o.f beam/rod elements to 
compute bending, torsional, axial, and in-plane deflection due to blade rotation and static 
aerodynamic loading. Test calculations on idealized model problems have been 
successful and correlation with measured deflections is awaiting full analysis of recently 
acquired data on the UH-60A model rotor hover tests. 

A combined linear programming/quadratic programming (LP/QP) approach to 
design optimization has been formulated and implemented. The optimization analysis 
accommodates all of the major planform variables as well as objective functions 
involving combinations of thrust, power (and its components), propulsive efficiency, or 
hover Figure of Merit User-imposed constraints can be imposed on these variables in 
addition to performance targets based on the selected objective functions. Sequential 
linear programming (SLP) is the basic mode of operation of the analysis, though an SQP 
option exists to find optimal solutions that exist away from the constraint boundaries. 
Limited updating of the tableau may be invoked to reduce CPU requirements in cases 
where the design evolution is well-behaved. 

A limited correlation study on helicopter rotors and tiltrotors in hover and axial 
flight has demonstrated that good accuracy can be achieved in integrated performance 
prediction using the EHPIC/HERO code. Demonstration calculations of the performance 
optimizer have also been carried out, using twist, chord, sweep, and anhedral in various 
combinations as design variables. In the cases examined, the optimum design obtained 
was a strong function of the constraints selected for that particular case (a well 
established feature of design optimization in general). The motivation for the constraint 
choices was to exercise the major features of the code in cases of practical interest; no 
claim is made as to the comprehensiveness of the tests undertaken here. However, two 
broad, interrelated trends have been identified that appear to constitute generally desirable 
design strategies. One of these is a manifestation of a generally well-understood feature 
of hover performance analysis: uniform distributions of bound circulation contribute to 
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the minimization of rotor power for a given thrust. A qualitatively similar principle 
appears to be at work in both the hover thrust maximization problems studied and the 
sample calculations in high speed axial flight. 

The second result is related to the first. Results from a wide variety of cases 
suggest that the addition of the 'down-up' twist or 'twist bucket' implemented on the UH- 
60A main rotor is a generally advantageous feature of blade design. The trend to this 
type of feature emerged in hovering rotor cases for combinations of twist with chord, 
sweep, and anhedral and for both low-twist, four-bladed rotors and high-twist, three- 
bladed tiltrotors. Preliminary calculations indicate this design feature evolves in power 
minimization problems even in cases where profile drag is excluded, suggesting that it is 
driven by the sensitivity of induced power to the wake-induced velocity field and the 
imperative to produce more uniform circulation distributions. Additional computational 
studies will be carried out to confirm this observation and to identify analogous strategies 
for other design variables. 

The results and computational experience of work to date on this analysis have 
also pointed out the need for additional effort in a variety of areas. The full exploitation 
of preliminary work on inclusion of models of vortex/tip loading using ANM would aid 
the resolution of vortex-induced loading. Further reductions in CPU time can be realized 
by exploiting existing simplifications available in the present optimization analysis. The 
optimized designs obtained with rigid blade analyses should be repeated with structural 
deflection in place to examine the effect on the solutions obtained to date. The 
preliminary work done on including constraints on the design from forward flight 
conditions should be extended to enhance the realism of the hover-optimized designs. 
Finally, extensions to include airfoil, blade thickness effects, and structural tailoring in 
the design optimization should be considered to supplement the already substantial 
capabilities of EHPIC/HERO. 
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APPENDIX A 


DESCRIPTION OF THE SIMPLEX ALGORITHM 

A brief summary of some of the features in the present implementation of the simplex 
algorithm is given here. Familiarity with the fundamental terminology and steps in an LP 
algorithm is presumed. For a review of the basic method see References 52 and 56. We 
shall concern ourselves here primarily with the extensions required for implicit 

incorporation of upper and lower bounds upon the state variables, AX, and for solving 
QP problems. In the following, we refer to those components of X which are not at one 
of their specified bounds as basic variables. The remaining entries form the non-basic 
set. 


The simplex algorithm employed here is a highly revised and extended version of the 
routine given in Reference 52. The routine seeks to solve the LP problem: 


Min J = J0 + c T AX 


(A- la) 


subject to the constraints: 

gj(X) + Vgj(X)AX * 0 

gk(X) + Vg k (X) AX = o 
AXj>0 

where 


fbj^AjjAX ,gj(x)<;o 
[b2^[A 2 ]AX , gj(X) > 0 
b 3 =[A 3 ] AX 


(A- lb) 

(A-lc) 
(A- Id) 


bl = -{ gj(X) 

and the matrices, 


h 2 = {gj(X) } ; h 3 = { g k (X) * sign{g k (X))) } 

(A-le) 


[A{1 = V gj (X) 


[A 2 ] = -V gj (X) ; [A 3 ] = -Vg k (X) *sign{Vg k (X)} 

(A- If) 


The ]2i are defined so as to have all entries ^ 0. The optimal solution is attained in two 

stages. In the first stage a feasible solution AX satisfying the imposed constraints is 
obtained by the introduction of slack variables and the construction of an auxiliary cost 

function. At the end of this stage, the actual optimization process is performed and AX 
varied via a sequence of pivot operations so that during each operation the cost is 
reduced. The process continues until no further reduction in cost is possible and thus 

returns with the optimal vector, AX- 


A.1 Slack Variables and Determination of a Feasible Solution 

The standard linear programming problem is posed in terms of a set of equality 
constraints and the requirement that all Axj^O. Thus the inequality constraints are 
transformed into equality constraints by the introduction of slack variables, $i and sj: 

i2i = [Aj] AX + fii ; b2 = [A2]AX-S2 (A-2a) 
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51 , 52^0 


(A- 2 b) 


Observe that the first constraint, bi = [A|] AX + Si, is easy to satisfy by setting, AX= 0 , 
and Si=hi. This is due to the positivity of the components of bl- This is not generally the 
case for the remaining inequalities associated with b2 and b3- Thus, in order to determine 
a feasible solution the constraint parameters, £ 2 and £3, are introduced to obtain the 
augmented set of constraints: 

bi = [Ail AX + Si ; b2 = [AJ AX - S2 + Z2 '■> 1*3 = [A3] AX + £3 (A-3a) 

Axj >0 ; Si.S 2^0 ; z 2 ,z 3 >0 (A- 3 b) 

Specification of an initial solution to Eqs. (A- 3 ) is now rendered trivial and is given by: 

ax=o ; si=bi ; 22=0 ; 53=° ; 22=^2 ; 23=^3 (A- 4 ) 

Since the bj^O. the and zj> 0 . Furthermore, {si, z 2 , £3} forms the set of basic (non- 
zero) variables and (AX, S2> S3} forms the set of non-basic variables. One now has a set 
of equality constraints, Eq. (A- 3 a), subject to conditions, Eq. (A- 3 b) which in conjunction 
with an appropriate linear cost function, define a linear programming problem. To obtain 
a feasible solution to the original problem, Eqs. (A-l), we must first remove the 
constraints variables, £i, i.e., conduct pivot operations aimed at reducing zi=0, or 
equivalently causing the £j to leave the basic set. Thus, the appropriate cost function to 
be minimiz ed is chosen as the sum of the components of £ 2 and £3: 


*z = £ z 2i + X z 3i (A-5) 

m2 m3 

where m 2 and m 3 are the dimensions of the vectors Z2 and I3 or equivalently the 
dimensions of b 2 and b3 respectively. The simplex machinery involving pivot selection 
and pivoting operations may now be applied to this minimization problem. If the 
minimizing solution subject to the constraints satisfies min{f z }=0 then this implies that 
the z 2 j and Z3} are zero and the corresponding original constraints are satisfied, i.e., b2 = 

[A-J AX - fi2» and b3 = [A3] AX- The simplex pivoting strategy guarantees that the other 
constraints involving bi and the positivity requirements upon the variables is maintained. 
If min{f z }> 0 , then one or more of the z 2 jor are positive and cannot be reduced further 
implying that the corresponding constraint is not satisfied. The initial value of f z 

f z = X h2i + X ^ ^ A ‘ 6 ^ 

m2 m3 

so that the progress of the simplex algorithm may be monitored by tracking the value of 
f z through the sequence of simplex iterations. The function, f z , must be expressed in 
terms of the non-basic variables which is easily accomplished by referring to the 
associated constraints: 

Z2=i22 - [AJAX + S2 ; 23= te-[A 3 ]AX (A-7) 


so that: 
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(A-8) 


/ \ 

f z = X b 2i + X b 3i“ X[ A 2] + X[ A 3] A 2£ + X s 2i 

m2 m3 \m2 m3 J m2 
This is conveniently summarized in the form of a tableau: 

Table A- 1 


Basic Variables 


£l 

%2 

23 


where, f z 0 = E b 2i + Xb3i (A-9) 

m2 m3 

and £„, denotes the vector of dimension m with unit entries, and [I m ] is the unit matrix of 
order m. The bottom row may be considered as an additional equality constraint as far as 
the simplex pivoting operations are concerned. The last three columns in the table 
correspond to the basic variables and form a unit matrix of dimension (ml+m2+m3). 
Since pivoting operations essentially permute the unit column vectors within the tableau 
it is only necessary to store the remaining tableau in computer memory and keep track of 
the basic and non-basic variables using index lists. Thus the last three columns are not 
actually stored in the array used in EHPIC/HERO, which leads to an m by m reduction in 
memory requirements with m being the total number of constraints. Furthermore, since 
the column vectors in the tableau for £2 are identical to those for 22 except for sign, and 
the pivoting operations affect both columns in the same way, an m by m 2 reduction in 
tableau space is effected by tracking these also with an index list, due consideration being 
made for sign. Procedurally, when a given component zjj leaves the basic set (i.e. 
becomes zero), it is replaced by the associated S 2 i variable, and the tableau modified by 
reversing the sign of the associated column vector. Thus, zjj is eliminated from the set of 
variables under consideration. When a component Z 3 i leaves the basic set it becomes a 
non-basic variable identified with column, j, in the tableau. Since Z 3 i is not used in 
subsequent computations, the column is henceforth omitted from consideration as a pivot 
column. 


A.2 Implementation of Upper and Lower Bounds 

In the above discussion a lower bound was implicitly imposed upon the Axj by the 

positivity requirements, Axj^O. This can easily be extended to any set of lower bounds, lj, 
by simple linear transformation of the defining variables. To impose upper bounds upon 

AX the pivoting strategy is modified following the procedures outlined by Luenberger, 
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Reference 56. Whereas in the original algorithm, the variables are referred to their lower 
(zero) bounds, each variable can now also be referred relative to its upper bound, 
UBND(i). An integer array, IBND(i) is employed to indicate which bound the variable Xi 
is referred to. The pivot selection procedure is now as follows. The upper bound for 

each variable, xj, is set by the user (design variables, d) or elsewhere in the code tec, 30- 
The upper bound for each of the fii, £ 2 , 22 , and Z3 is arbitrarily set to a large number. The 
same is true of the Lagrange multipliers when conducting QP problems. 

1) The non-basic variable to be pivoted, Xj, is determined. This is done in the same 
manner as in the single bounded case and is accomplished by searching the row 
corresponding to f z (when finding a feasible solution) or J (when determining the optimal 
solution) and finding the column, j, with minimum entry. If this value is positive then no 
further improvement is possible. Otherwise we proceed with the step 2). 

2) There are now three possible moves a) The non-basic variable, Xj, simply goes to 
its opposite bound; b) xj enters the basic set by being pivoted with a basic variable, xj, 
and xj returns to its old bound (this is the usual pivoting procedure for the single bound 
case); c) same as b) except that the variable Xj is brought to its opposite bound. To 
determine both the action taken and for cases b) and c) which basic variable, xj, is to be 
pivoted determine the minimum of: 


a) UBND(i) 

b) min a i0 /aij 
a ;i >0 


c) min ( ajo -UBND(i) )/ay 

a ij < 0 


where ay is the entry in the above table corresponding to constraint, i, and nonbasic 
variable, j, and aio is the leftmost column entry of constraint, i. 

3) Depending upon the minimum value in step 2), execute the operations associated 
with a), b), or c). For details of the pivot operation see Reference 56 noting that this is 
implemented upon the compact storage scheme of Reference 58. 


A.3 Extension to QP Problems 

Each of the quadratic programming problems is solved by a modified version of the 
simplex algorithm. The Kuhn-Tucker conditions are simply linear equality conditions 
and if these are satisfied together with the original imposed constraints then one has the 
solution to the QP problem. However, the difficulty lies in meeting the criterion, Eq. (5- 
14b,c), which forms a nonlinear constraint The solution approach follows that of Wolfe 
(Reference 59) and involves a two stage procedure as in the LP case. In the first stage a 
feasible solution satisfying the original constraints is derived in an entirely analogous 
manner as in LP. Starting with the feasible solution thus obtained, the additional Kuhn- 
Tucker equality constraints are then imposed and a new feasible solution satisfying the 
entire set of constraints is sought The nonlinear condition, Eq. (5-14b,c), is implicitly 
satisfied by modifying the logic used in selecting the pivot element in the simplex 
pivoting process. 

The Kuhn-Tucker criteria augment the original set of imposed constraints by n further 
equality constraints together with the orthogonality requirements upon pairs of Lagrange 


114 



multipliers and constraints. The equality constraints fall directly into the LP capability of 
the simplex algorithm. It is the enforcement of orthogonality between certain pairs of 
variables necessitates modifications of the method. Recalling the Kuhn-Tucker 
conditions, Eqs. (5-14): 


Q = VJ(X) + [B] AX+[aJ' -Aj]ij+[l„](ii 2 -/£,)+[A3]it (A-10a) 


Mli (Axj Alj) — 0 , fi2\ (Auj Ax|) — 0 , /tjj , fJ.2i ^ ® » (1 l*****n) 

(A- 10b) 


AjSj = 0 ; Aj £ 0 (j=l,...,mi) 

Ak unbounded (k=mi+l,...,me) 


(A- 10c) 


(A-lOd) 


where {sj}={si , £2} are the slack variables associated with the inequality constraints, and 
JaJJ = [V g k (X)] T and is equivalent to A3 T except for the sign reordering defined by 
Eq. (A- lb). There are n new equality constraints and m+2n new Lagrange multiplier 
variables. The condition Aj Sj = 0 is equivalent to the stipulation that Aj can only enter 

basic set if Sj is zero implying that the corresponding inequality constraint is active (i.e., 
satisfied exactly). 

The Lagrange multipliers associated with the imposed equality constraints is 
unbounded and so can be eliminated from the set of variables by solving for them directly 


using m3 of the n Kuhn-Tucker equality constraints. The matrix. 


[Aj] = 


A l 3 ' 

— T 

An-m3. 


is partitioned, 


(A- 11) 


where the rank of [a£ 3 ] is m3. The rows of [aJ 3 ] are chosen by searching for 
maximum pivots so that [a^ j is well conditioned. Then, 

ik =-[A m3 r T {vj(X) + [B]AX + [Af -Aj]ij + (/£ 2 -M 1 )} n]3 

(A-12) 

where the subscript (*) m 3 refers to the row partitions of the argument corresponding to 

those of [aJ 3 1. These equations may now be dropped from the Kuhn-Tucker 
conditions. The remaining n-m3 Kuhn-Tucker equations are modified by substituting for 
Ay . The reduced set of n-m3 Kuhn-Tucker conditions then takes the form: 


Q = -£+[Q]AX + [R]Aj + T(u 2 -Ui) (A- 13a) 

ci>0 (A- 13b) 
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together with the original orthogonality requirements, Eqs. (A-10b,c). The matrices [Q], 
[R], and [T] are obtained by substituting for in the remaining n-m3 rows of the original 
Kuhn-Tucker equations, Eq. (A- 10a). The condition q £ 0 is effected by reversing the 
sign of each entry in row i if necessary. 

In the first stage, the additional vectors, 12 and 23. are temporarily introduced as 
before so that an initial feasible vector satisfying all augmented imposed constraints, 
Eq. (A-3), can be found. Similarly, the vector, zkt is introduced so that the Kuhn-Tucker 
equality constraints become: 

£ = [Q]AX + [RU. + T(u 2 -ui) + ZKT (A-14a) 

acr^Q < A -' 4b) 

A feasible solution vector satisfying all constraints and the orthogonality requirements is 
obtained by simply setting all of the Lagrange multipliers to zero, zkt=£. and initializing 
the remaining variables as in Eqs. (A-4). The resulting tableau reads: 


Table A-2 



[AiJi H 

[0] 

mm 

mm 

[0] 

Si 


r [Aji 

- [Wl 

[0] 

[0] 

[0] 

JSL 

B9 

[A 5 Ji 

[0] 

mm 

HI 

[0] 

*3 




mm 

[RJ 

msm 

mum 

-KT 

^z0 

m2 m3 

S( A 2]ki + X[ A 3]ki 

k k 

-{Cm2) 






where the final row is formed in exactly the same manner as in the LP case. The last 
three blocks of the final row are not used in stage 1. During this first stage, the Lagrange 
multipliers remain at their lower bounds and are skipped when determining the pivot 
columns. Therefore, the orthogonality requirements imposed by the Kuhn-Tucker 
conditions is maintained. Note however, that the pivot operations basically add multiples 
of the pivot row to all of the other rows so that the matrix [Q] will be modified. The 
matrices [R] and [T] remain unaffected since all other entries in the associated column 

blocks are zero. 

From the above table it is clear that the columns pertaining to Pu and p 2 j identical 

except for the sign. Also, it is clear that it is impossible to have both Pu and p 2 i both 
belong in the basic set since that would imply that the associated xj variable is 
simultaneously at its upper and lower bound. The row manipulations taking place in a 
pivoting operation modify the columns associated pjj affect pjj in identical manners. 
Thus it is not necessary to represent both variables in memory since knowledge of one of 
the columns implies the same for the other. In practice, an integer array is used to 

ascertain which of the variables, Pn or p 2 i is currently represented in the tableau. When 
a column pertaining to a Pu or p 2i variable is considered during the pivot selection tests, 
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it is e xamin ed twice, first with the unaltered column and then again with all entries in the 
column reversed in sign. Thus an n by n reduction in memory space is realized. 

The second stage tackles the Kuhn-Tucker equality constraints. Specifically, a 
feasible solution satisfying these equality constraints is sought. In an analogous fashion 
to the construction of f z in Eq. (A- 8 ) which was used to obtain a feasible solution for the 
imposed constraints, an auxiliary linear cost function is created, 

n-m3 

fKT= E(acr)j (A-15) 

i=l 

As before, this is expressed in terms of the non-basic variables by substituting for the 
ZKT- If one examines the tableau, Table A-l, and the expression for f z in Eq. (A- 8 ), it is 
seen that the constant term, f z o, is simply the column sum of the left-most column with 
the summation range extending over those rows corresponding to the Z 2 i and Z 3 i. 
Similarly, the factor associated with a given non-basic variable, xj, in the expression for 
f z can be obtained simply by summing the entries in the column above it over the same 
summation range. Or, 


ml+m2+m3 

( a z)j= X a ij (j=0, n+ml) 
i=ml+l 


(A- 16) 


where fa is the row in the tableau associated with f z , j =0 corresponds to the leftmost 
column, i=0 corresponds the row associated with J, and j=n+ml corresponds to the 
righmost non-basic variable, in this case (si) m l- The same procedure is carried out for 
the fKT with the row entries being computed as a simple sum of the associated column 
entries. The range of summation now extends over the n-m3 rows identified with the 
transformed Kuhn-Tucker conditions (i.e., subsequent to the elimination of m3 rows 

when solving for the unbounded 2dc). Furthermore, the column range, j, is extended to 
n+ml+n+mi so that these sums are also formed for the remaining n+mj Lagrange 
multipliers which at this point arc still non-basic variables set to their initial zero values. 


fz 


S b 2i + S b 3i“ 
m2 m3 


E[A2]+L[A 3 ] 


Vm2 


m3 


AX + Xs 2i 

m2 


At the end of the first stage, the tableau has the form: 


(A- 17) 


Table A-3 
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I <4 
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mm 



k=m+l 

k=m+l 

HU 

■ 



where (•)* denotes that the entries have been altered due to row manipulations during 
pivoting. Here, the scalar, 

m+n-m3 

f KT0 = X c k (A* 18 ) 

k=m+l 

The second stage then is concerned with minimizing fj^r- If the optimal value of fKT=0> 
then by implication all the zkt= 0 thus verifying that a feasible solution to the complete 
problem of imposed constraints plus the Kuhn-Tucker conditions, or equivalently, the 
optimizing vector of the QP problem, has been found. This presumes that the 
orthogonality conditions, Eqs. (A-10b,c), are maintained throughout, and this must be 
explicitly enforced in the second phase. To this end, a battery of logical tests is 
conducted during selection of the pivot column. For each candidate column, the 
associated non-basic variable together with its set of conjugate variables (i.e. the second 
in the pair of variables that define the orthogonality requirement) is found. For any xj the 

set of conjugate variables is |iij and |i2j; for slack variable, sij or S2j, the conjugate 

variable is the corresponding entry of 2tj; etc. The bookkeeping is somewhat involved 
due to the compact storage of the simplex tableau and the implicit incorporation of lower 
and upper bounds upon the variables, but in essence the associated conjugate variables 
are tested for membership in the basic set. If the membership structure is found to be 
incompatible with the orthogonality stipulations, Eqs. (A-10b,c), then that column is 
skipped. If no column can be found that both satisfies the orthogonality conditions and 
decreases the auxiliary cost, fjcr* then the possibility of a non-basic variable pivoting with 
one of its own conjugate variables is examined. For example, one could pivot the non- 

basic Xj component with its corresponding Jijj assuming that p.ij presently belongs to the 
basic set. Such situations arise rarely in practice, and pivot columns satisfying the 
orthogonality conditions are usually available. 
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